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ABSTRACT 



A local geoid model to predict the geoid heights in the vicinity of Monterey 
Bay, California, was developed to use Global Positioning System (GPS) differential 
positions and known Mean Sea Level (MSL) with the method of collocation. The 
local geoid models were based on Rapp's 360 degree x 360 order global geoid 
model determined from gravity measurements. Control data were adjusted by least 
squares to solve for the parameters in the local geoid model. Also studied were 
factors that affected the GPS-measured ellipsoid height differences. These included 
(1) comparing GPS differencing solutions, (2) standard error of GPS observations, 
(3) corrections for surface meteorological values, and (4) observation durations for 
GPS. 

The data used in this research were taken from GPS measurements on the 
campus of the Naval Postgraduate School (NPS), an area about 100 m x 630 m and 
in an area approximately 15 km x 33 km near Monterey, California. The time 
period was from February 5, 1988, to May 12, 1988. 

The accuracy of the predicted geoid heights is ± 2 cm if a six-parameter model 
is used for the large area, and ± 2 to 10 mm if a five-parameter model is used for the 
NPS campus. 
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I. INTRODUCTION 



The Global Positioning System (GPS) is able to establish precise relative 
positions in the World Geodetic System of 1984 (WGS 84). A Trimble GPS 
receiver, which has the capability of measuring carrier phase, was used to 
determine the vector base line in space. The components of the base line are 
expressed in terms of cartesian coordinate differences (AX, AY, AZ) [Remondi, 
1984]. These vector base lines can be converted to distances, azimuths and the 
ellipsoid height differences, Ah, relative to the WGS 84 Ellipsoid. 

The results of several tests and operations have clearly shown that GPS 
survey methods can replace conventional horizontal survey methods. 
Comparable accuracies have also been achieved for GPS-derived ellipsoid height 
differences, Ah. The problem of converting these ellipsoid height differences. 
Ah, to orthometric height differences, AH, remains to be resolved. For example, 
in engineering surveying applications WGS 84 coordinates must be transformed 
to the North American Datum of 1983 (NAD 83) system. The GPS obtains 
ellipsoid heights, rather than the orthometric heights; the geoid height, N, must 
be calculated to obtain the latter. One of the problems in this transformation is 
the accurate determination of the local geoid height, N. 

For GPS survey applications, a geoid model should provide geoid heights 
with an accuracy commensurate with that of the ellipsoid height, h, so the 
accuracy of the derived orthometric heights is not reduced. In the future 
differential positioning will be widely used in GPS surveying, and, therefore, 
only the geoid height differences between stations will be required. 

Geoid height computation techniques include the following: 

(1) Geoid height differences in the U.S. can be determined from gravity 
data and the Stokes’ integral method, or from astrogravimetric data and least 
squares collocation methods. These methods lead to uncertainties that are 
typically 1 to 10 cm for distances of less than 20 km and 5 to 20 cm for distances 
between 20 to 50 km [Zilkoski, 1988] 

(2) A comparison of a data set from GPS with gravimetrically determined 
geoid heights using least squares collocation techniques shows discrepancies 
between the two data sets of about ± 2 cm for a maximum intersection distance of 
approximately 50 km [Denker and Wenzel, 1987]. 
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(3) Mean gravity anomalies, deflections of the vertical and a geopotential 
model calculated to degree and order 180 have been used to determine geoid 
heights in the area bounded by (34° < (p < 42°, 18° < X < 28°) [Tziavos, 1987]. 
The method used was that of least squares collocation. By using empirical 
covariance functions for the data, suitable weighting functions for the different 
sources of observations, and the optimum cap radius around each point of 
elevation, an accuracy better than ± 0.60 m was obtained for geoid heights. 

The main objective of this thesis is to use a model that predicts N in a region 
near Monterey, California, from the GPS differential positions and known mean 
sea level (MSL) using the method of collocation. Also studied were factors that 
affected the ellipsoid height differences, Ah, obtained from GPS measurements, 
such as comparing GPS differencing solutions, standard error of GPS 
observations, corrections for surface meteorological value , and observation 
durations for GPS measurements. The data used in this research were taken 
from GPS measurements on the campus of the Naval Postgraduate School (NPS) 
an area about 100 m x 630 m and in an area approximately 15 km x 33 km near 
Monterey, California. The time period was from February 5, 1988, to May 12, 
1988. 

The accuracy of the predicted geoid height is ± 2 cm if a six-parameter 
model is used for the large area, and ± 2 to 10 mm if a five-parameter model is 
used for the NPS campus. 
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II. GEOID HEIGHT 



A. THE GEOID 

Surveyors and engineers, in most cases, are interested in the orthometric 
height, H, as measured above the reference surface of the geoid. The ocean is 
considered to be freely moving, homogeneous and only subject to the force of 
gravity. When a state of equilibrium is achieved, the surface of this idealized 
ocean assumes a level surface of the gravity field. It may be regarded as also 
extending under the continents. This level surface is called the geoid. If the 
potential, W, is given as a function of the coordinates r, <{> and X, then the geoid is 
given by [Moritz, 1984] 

W = W(r,<j>,A.) = W 0 

The geoid is a closed and continuous level surface which extends partially 
inside the solid body of the earth. The direction of the gravity vector at any 
point (plumb line or vertical) is normal to the geoid. The curvature of the geoid 
displays discontinuities at abrupt density variations. Consequently, the geoid is 
not an analytic surface, and therefore not a practical reference surface for 
position determinations. The geoid however, is well suited as a reference surface 
for potential or height differences, which are obtained by precise levelling in 
combination with gravity measurements. 

To establish the geoid as a reference surface for heights, the ocean water 
level is recorded and averaged over long intervals (> 1 year) using tide gauges. 
The MSL thus obtained represents an approximation to the geoid. The National 
Geodetic Vertical Datum of 1929 (NGVD 29) was derived for land surveys from 
a general adjustment of the first order levelling net of both United States of 
America and Canada. In the adjustment MSL was observed at twenty-one tide 
stations in the United States and five in Canada. The geoid established by this 
method may deviate by ± 1 to ± 2 m from a level surface due to periodic, 
nonperiodic and secular variations [Torge, 1980]. 
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B. THE WGS 84 ELLIPSOID 

The development of WGS 84 [DMA, 1987] was initiated by the United States 
Department of Defense for navigation and weapon systems. The Defense 
Mapping Agency (DMA) developed WGS 84 as a replacement for WGS 72. The 
defining parameters and reference frame orientation of the WGS 84 Ellipsoid 
and the WGS 84 Ellipsoid Gravity Formula are those of the internationally 
sanctioned Geodetic Reference System of 1980 (GRS 80). Accordingly, a 
geocentric equipotential ellipsoid is defined by the semimajor axis (a), the 
flattening (f), the equatorial gravity (y a ), and the angular velocity (co). The WGS 

84 Ellipsoid used the values: 



The reference system for GPS is WGS 84. The precise geocentric 
coordinates obtained from GPS receivers are in WGS 84. 



C. THE ORTHOMETRIC HEIGHT, H, THE ELLIPSOID 
HEIGHT, h AND THE GEOID HEIGHT, N 



In geodetic applications three different surfaces or earth figures are 
normally involved. First is the earth's actual topography; second, the geometric 
surface, or ellipsoid and; third, the equipotential surface, the geoid. The 
relationship between the earth's topography, the ellipsoid and the geoid in a 
section through the earth's surface is shown in Figure 1. 



f 

Ya 



a 



(0 



6378137 m 
1/ 298.257223563 
987.03267714 gals 
7.2921 15xl0' 5 rad/s 
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Figure 1. Relationship between the earth's surface, the geoid and 
the ellipsoid. 

Two features in the figure are of particular interest : 

1. The deflection of the vertical, 0, defined by Pizzetti as the angle at the 

geoid between the direction of the plumb line (normal to geoid) and 
the normal to the ellipsoid through the point, P, on the geoid [Torge, 
1980]. 

2. The vertical separation, N, between the geoid and the ellipsoid. 

The deflections of the vertical and the geoid heights, N, depend on the 
ellipsoidal coordinates and, hence, on the parameters of the reference ellipsoid 
and its position with respect to the earth. If they are referred to the 
geocentrically situated mean earth ellipsoid, then they are referred to as absolute 
quantities; otherwise, they are relative quantities. The absolute deflections of the 
vertical in flat terrain and the highlands assume values between one and ten 
seconds of arc; in mountainous areas, they vary between one-half and one minute 
of arc. As a result of density variations within the earth, the geopotential 
surfaces, including the geoid, have irregular shapes. Absolute geoid heights, 
however, rarely exceed 100 m [Torge, 1980]. 

The orthometric heights, H, shown in Figure 2 are referred to an 
equipotential surface, the geoid. The orthometric height of a point on the earth's 
surface is the distance from the reference surface to the point, measured along 
the plumb line normal to the geoid. The ellipsoid height, h, of a point is the 
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distance from the reference ellipsoid to the point, measured along the line which 
is normal to the ellipsoid. For purposes of simplicity, the orthometric height, H, 
and the ellipsoid height, h, are shown to be along a common vertical. In most 
cases, this would cause a very small error that is considered insignificant 
compared to present uncertainties of the geoid height N estimates. The geoid 
height, N, is defined: 



N = h - H 




Figure 2. Relationship between the geoid height, N, the ellipsoid 
height, h and the orthometric height, H. 

The orthometric height has greater physical meaning than the geometrical 
ellipsoid height. The orthometric height has traditionally been determined by the 
technique of levelling in which increments of height are obtained from the 
intersection of the line of sight of a level instrument, tangential to the 
geopotential surface and passing through the level axis, with two graduated rods. 
Accurate orthometric height information is needed for precise engineering 
operations such as the construction of dams, pipelines, and tunnels. 

Geoid heights have been accurately determined for some major geodetic 
datums such as the North American, European and Australian. These datums are 
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well supplied with astrogeodetic deflections and have fair gravity coverage. The 
standard error of relative geoid heights in these areas is about 2 or 3 m. 



D. THE GEOID HEIGHT FROM GRAVITY MEASUREMENT 

Both the U.S. National Geodetic Survey (NGS) of National Oceanic and 
Atmospheric Administration Service (NOAA) and Trimble Navigation versions 
of Rapp’s 360 degree x 360 order model (OSU86F) were available to me. Rapp 
computed two potential coefficient fields that are complete to degree and order 
360 [Rapp, 1986]. One field (OSU86E) excludes geophysically predicted 
anomalies, while Rapp’s other model (OSU86F) includes such anomalies. These 
fields were computed using a set of 30-minute mean gravity anomalies derived 
from satellite altimetry in the ocean areas and on land from standard 
measurements. 

Gravity anomalies can be observed and then used to compute the geometric 
deviation of the geoid from the ellipsoid. The expression for the gravitational 
potential is written in the following form [Rapp, 1986]: 



... kM ~ 
V(r,(j),X)= — — 1 + X 
T L 1=2 



a 

rj 



X ( C, cosmX + S, sinmX, 1 P, 
m =o v lm lm ) lm 



(sin<}>) 



where 

r, <{>, X : geocentric coordinate 

kM : geocentric gravitational constant 
a : equatorial radius of the reference ellipsoid 

C , S : fully normalized potential coefficients 

P lm : fully normalized Legendre function of degree 1 and m 
y : normal gravity 

The potential at a point, U, is the scalar sum of V and centrifugal force potential 
V' [Ewing, 1976]: 



U(r,4>A) = V(r,4>,X) + V(r,4a) 



The difference between observed gravity potential W and the computed normal 
gravity potential U is denoted by T [Moritz, 1984], so that 
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W(r,<|>,X) = Uti,^) + T(r,<t>,X) 
compared the geoid defined by the potential W is given by 

W(r,<M.) = W 0 

A reference ellipsoid with the same potential, W = U , is given by 
U(r,<ta) = W 0 

A point P on the geoid is projected onto the point Q on the ellipsoid along the 
normal to the ellipsoid PQ = N. PQ is the distance between geoid and ellipsoid at 
the point and is called the geoid height (Figure 3). 




GEOID 
W = w 0 



ELLIPSOID 

u = w 0 



Figure 3. Geoid and ellipsoid (Moritz [1984, Fig. 2-12]). 



The disturbing potential T can be written as [Rapp, 1986]: 



tM °° 

T(r,<a) = v- I 
Y 1=2 



' X X Ct ¥,“(<]>, X) 

m=0 a=0 
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where 



C^= {C lm , a = 0 , and S Jm , a = 1 } 

Y^(<{),^) = {P^sin^cosm^, a = 0, and P^sin^sinmT., a = 1 } 

Since the N is relatively small compared to the ellipsoid geocentric radius R', 
then [Ewing, 1976] 



and 



T = yN 




where y is the normal gravity at (<{), X). 

The gravity anomaly, Ag, in the Molodensky surface free-air anomaly sense 
is given by [Moritz, 1984]: 

Ag p ~ §p ' Yq 



The boundary condition that relates Ag and T is [Moritz, 1984] 



Ag 



0T l dT 
3h + r 9h 



where h is the distance along the plumb line direction. Neglecting deflections of 
the vertical we have [Rapp, 1986] 



Ag 




oo 

1 ( 1 - 1 ) 

1=2 




i £< cL-eK-ea-es)^) 

m=0 a= o 



where ^ ( i = h, y) are ellipsoidal corrections. 

The Stokes function, S(\|/), can then be used to solve the geoid heights above the 
geodetic ellipsoid [Ewing, 1976]. 
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N g = 4n~y Jo "Jo Ag(V,a)S(\|/)sin\|/d\|/doc 

'm 

where 

Y m : the mean value of normal gravity 

\j f : the angular distance between the point where N is being 

determined and the area where the effect of Ag is being 
considered 

a : the azimuth from the affected point to that causing the effect 
Ag : the gravity anomaly 
S(y) : the Stokes function 

The Stokes function is given by [Ewing, 1976] 

S(\j /) = esc + 1-6 sin ^ - 5 cos\j/ - 3 cos\|/ ln(sin + sin 2 

The gravitational potential using degree 360 and order 360 has been tested 
through comparison of Doppler station geoid heights with geoid heights from 
Rapp's versions of geopotential models. The agreement between the two geoid 
height measurements is approximately ± 1.6 m [Rapp, 1986], 
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III. GEOID HEIGHT FROM GPS 



A. THE GLOBAL POSITIONING SYSTEM 

The Global Positioning System (GPS) calls for a precise navigation system 
divided into three segments: space segment, control segment and user equipment. 
The space segment will consist of three orbital planes of satellites at inclinations 
of 120° in circular orbits at altitudes of 20,000 km. Each plane will eventually 
contain six to eight satellites to give the three dimensions of position, velocity and 
precise time 24 hours a day anywhere in the world. The control segment consists 
of the ground stations necessary to track the satellites, monitor the system 
operation, and periodically provide corrections to the navigation and time 
signals. Each satellite broadcasts signals containing information on its position. 
The GPS satellite transmits signals at two L-band frequencies (1227 and 1575 
MHZ) to permit corrections for ionospheric corrections. The signals are 
modulated with two codes: P, which provides for precise measurement, and C/A, 
which permits easy lock-on to the desired signal. The user segment consists of 
the equipment necessary to convert the satellite navigation message into useful 
navigation information. 

The navigation message contains the data that the user's receiver requires to 
perform the operations and computations for successful navigation with GPS. 
The data includes: (1) information on the status of the Space Vehicle (SV); (2) 
time synchronization information for the conversion of the C/A to P code; (3) 
parameters for computing the clock correction and the ephemeris of the SV; and 
(4) corrections for delays in the propagation of the signal through the 
atmosphere. In addition the data contain almanac information to define the 
approximate ephemeris and to give the status of all the other SV information 
which is required for use in signal acquisition [Milliken, 1980]. Ranges to the 
satellites are determined by signal transit times multiplied by the speed of light 
(299,792,458 m/sec). The transmitted message contains ephemeris parameters 
that enable the user to calculate the position of each satellite at the time of the 
transmission of the signal. 
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B. CARRIER PHASE MEASUREMENTS 

GPS measurements can be made using the pseudo-range and the carrier 
phase. The pseudo-range is essentially a measurement of distance contaminated 
by clock error. When four satellites are observed simultaneously, the three 
dimensional position of the ground receiver can be determined along with the 
receiver clock offset at a single epoch. The accuracy of pseudo-ranges is affected 
by multipath effects, which depend on the antenna design, and its height above the 
ground. 

Carrier phase measurements are more precise than pseudo-range and are 
not as vulnerable to multipath effects. They can be used to compute the precise 
base line components AX, AY, AZ between two receivers. Phase measurements 
are made by beating the received carrier with the signal from a local oscillator 
internal to the GPS receiver. The slant range from a GPS receiver to a satellite 
can be modelled in terms of time. It takes the signal time to travel between the 
satellite and the receiver or the equivalent number of cycles. The range of the 
cycles will consist of an integer and fractional number of cycles. When a 
receiver locks onto the carrier signal, it can immediately measure the fractional 
part and begin counting subsequent integer cycles, but it can not measure or 
account for the initial integer number of cycles that preceded the initial fractional 
part. The initial integer ambiguity which biases the subsequent measurements is 
called the initial integer ambiguity bias. 

GPS uses a one-way carrier beat phase. The GPS satellite and receiver are 
controlled by separate clocks. The satellite clock generates the signal, and the 
receiver clock detects when the signal arrives. An error in the synchronization 
of the clocks of 1 microsecond creates an error in range of 300 m. 

C. ONE-WAY CARRIER PHASE MEASUREMENT 
DIFFERENCING 

A single difference, SD(j,i), is formed by differencing carrier beat phase 
observables from two receivers 1 and 2 at the same observation epochs i of same 
satellite j. The equation is given by [Remondi, 1984]: 

SD(j,i) = S(2,j,i) - S(l,j,i) 
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where S(k,j,i) is the raw, unpreprocessed, fractional phase plus the count made at 
epoch i by receiver k for satellite j. The main advantage of the single difference 
is that it reduces or eliminates satellite orbital and clock errors, because they are 
common to both receivers. Its disadvantage is that one can not exploit the integer 
nature of integer ambiguities. Thus, for short base lines, the ultimate in accuracy 
may not be achievable [Remondi,1985], 

A double difference, DD(j,k,i), is formed by differencing single differences 
between a reference satellite j and another satellite k at the same epoch i. The 
equation is given by [Remondi, 1984]: 

DD(j,k,i) = SD(k,i) - SD(j,i) 

The advantage of the double differences is that the receiver clock dependent 
terms are eliminated because the differences for each epoch are correlated. The 
significance of the removal of those terms is to reduce from nanoseconds to 
microseconds the timing accuracy required to achieve one cycle accuracy. For 
short base lines, the integer ambiguities can be isolated, since the contribution 
made by the clock drift is reduced [Remondi, 1985]. The Trimble 4000SX 
receiver achieves sub-microsecond accuracy by using the C/A code timing 
information [Ashjaee, 1985]. 

A triple difference, TD(j,k,i), is formed by differencing the double 
differences for the same satellite pair at some integer of succeeding epochs i+1 
[Remondi, 1984]. 

TD(j,k,i) = DD(j,k,i+l) - DD(j,k,i) 

The advantage of the triple difference is that it eliminates all the time independent 
terms, namely the initial integer ambiguities, and becomes insensitive to the 
initial ambiguities and any cycle slips when, the receiver loses lock. The 
disadvantage of the triple difference is, another level of correlation, loss of 
resolution and a greater number of observations. Triple differences are already 
correlated with respect to satellite because of the underlying double differences 
and are further correlated with respect to time because consecutive triple 
observations will have common DD(j,k,i+l) terms. 
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For short base lines local area integer ambiguities can easily be resolved 
because unmodelled errors are highly correlated between the two antenna sites and 
are mostly eliminated by differencing. Algorithms can take advantage of the 
integer nature of the initial ambiguities and solve for them [Remondi, 1984]. 

D. GEOID HEIGHT FROM GPS 

Surface fitting techniques can be used with GPS-derived geoid heights. GPS 
stations are likely to be close together, of the order of a few tens of km, and the local 
geoid can be estimated directly if levelling data is available in the area. This 
together with it's relative simplicity makes the method practical for correcting GPS 
heights. 

It is possible to use the two sets of elevations (that is, levelled and GPS- 
derived) to define two distinct planes. The published levelled elevations are 
referred to the geoid and the GPS elevations are referred to an ellipsoidal surface. 
If both elevations are made equal at one bench mark (control station), then in 
general, the other bench marks will have two elevation values. After several 
different models were studied two mathematical surface models were chosen in this 
study. The five-parameter model found suitable for use on the NPS campus is 

N 0 (X,Y,Z) = h 0 + Ah - H + aAY + bAX 2 + cAY 2 + dAXAY 

and the six-parameter model selected for a larger area near Monterey, California, is 
given by 



N 0 (X,Y,Z) = h 0 + Ah - H + a' AX + b'AZ + c’AX 2 + d'AY 2 + e’AXAY 



where 

N 0 (X,Y,Z) : the global geoid height, obtained by gravity measurements 
h 0 : the ellipsoid height of the reference point, including a constant 
correction to N 0 (X,Y,Z) at the. reference station 
Ah : the ellipsoid height difference with respect to the reference station 
H : the published or levelled elevations at each control station 
a, b, c, d, a’, b', c', d\ e' : the coefficients to be determined 
AX, AY, AZ : the coordinate differences in WGS 84 
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These models can be solved to determine the east- west and north- south tilts 
which are absorbed by two rotations, one around the north axis (AY) and the 
other around the east axis (AX) in the horizon system. Their separation is 
absorbed by the scale correction [Zilkoski, 1988]. The local geoid heights, N, 
can be found from the equation 

N = N 0 (X,Y,Z) + AN 

where AN is the variation of geoid height in local area. The equation is given by 

AN = H + N 0 (X,Y,Z) - Ah 

AN = h 0 + aAY + bAX 2 + cAY 2 + dAXAY 
or 

AN = h 0 + a'AX + b'AZ + c'AX 2 + d'AY 2 + e'AXAY 

To solve for h 0 , a, b, c, d, a', b', c', d' and e', the global geoid heights 
N 0 (X,Y,Z) can be obtained from the Geoid.exe program described in Chapter 
IV. The geoid model can be rearranged to give 

H = h 0 + Ah - N 0 (X,Y,Z) + aAY + bAX 2 + cAY 2 + dAXAY 
or 

H = h 0 + Ah - N 0 (X,Y,Z) + a’AX + b’AZ + c’AX 2 + d’AY 2 + e’AXAY 

Then h 0 a, b, c, d, a', b', c', d' and e' can be solved by the least squares method. 
The geoid height, N, found by this method appears to be adequate for areas up to 
50 km x 50 km where the geoid is smooth [King, 1985]. 

E. METHOD OF COLLOCATION 

The method of collocation was derived from least squares interpolation. 
This method was used to predict the geoid height, N, for a local area. The geoid 
model is given by 

N = N 0 (X,Y,Z) + AN + S + n 
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where 



S : the signal 
n : the noise 

N 0 (X,Y,Z) + AN : the system function from the surface models at control 
points where both h and H are known 

Then, for a given N 0 (X,Y,Z) the method of collocation can be used to determine S q 
at these control points and to predict S p at other points in the local area. At any point 
in the local area the value of h can be determined by using differential GPS 
measurements between a known point and any other point by the equation 

h = h 0 + Ah 

The value of H at any point is given by 

H = h - N 0 (X,Y,Z) - AN - S - n 

The value of H at any point in the local area, which depends on the value of h at 
the control points and the differential GPS measurements, can be predicted with an 
accuracy of ± n. The general form of the observation equation in the method of 
collocation is [Jeyapalan, 1977]: 

x = A»X + S + n + 0*S 
q q p 

where 

x : the vector of the observation (x = Ah - N 0 (X,Y,Z) - H) 

A : a given rectangular coefficient 

X : the vector of the systematic parameters (h 0 , a, b, c, d, or h 0 , a', b', c', 
d\ e') 

S q : a signal vector at q observation points 
n q : a vector of measuring errors, noise at q points 
S p : a signal vector at p unknown stations 
• : indicates matrix multiplication 

0 : the null matrix 
If 

Z — S + n 
q q q 
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Then 



x = A*X + Z + 0*S 
q q p 

then 

X = (A T Cq~ 1 A )' 1 A T Cq' 1 x 
and 

S = C C’Vx-AX) 
p pq q 

where the variance-covariance matrix C of the Z and S vectors is 

q p 




The essence of this method is that by some means a covariance matrix can 
be assigned to the signal. For noise it will be possible to assign a diagonal weight 
matrix. 

S p are the values of the signal at the interpolated stations. Suppose there are 
q observations (and values of S q ), p interpolated values of S p and m model 
parameters. The covariance matrices are the following [Bomford, 1980]: 

(i) C q , the expected covariance between the observed x's for all pairs of 

the q observations. It is a q x q matrix. 

(ii) C sq , the expected covariance between the signals for all pairs of the q 

observations. It is a q x q matrix. 

(iii) Cn q , the expected value of n at each station. It is a q x q diagonal 

matrix. 

(iv) C pq , the expected covariance between all pairs of mixed observed 

and interpolated signals. It is a p x q matrix. 

(v) C , the expected covariance between the signals at pairs of 

interpolated stations. It is a p x p matrix. 

The variance of the noise, C nq can be estimated in the usual way according to 
the circumstances at each station, different types of instrument, etc. The noise at 
different points is independent; hence 
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nq 
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where G n is the standard error of the noise. 

The expected covariance between the observation x’s, C q can be obtained 
by 



C = C + C 

q nq sq 

because the signal is small. 

The C sq can be be computed from a simple function whose parameters can 
be determined by using empirical data. The covariance function is positive and 
definite. These two characteristics are found in many functions. Functions 
commonly used for covariance may be constant, sinusoidal, Gaussian, 
exponential, exponential cosine, and exponential sine and cosine. In this thesis the 
sinusoidal function was used! The function is given by 

C = B sin(kr) 

sq 

C(r) = C nq + B sin(kr) 

where C nq is the standard deviation of the control stations and B, k are 
coefficients to be determined. 

The parameters can be determined by least squares, and residuals at each 
point can be computed. From the residual, the covariance between points at a 
distance (r) can be computed by 



C(r) = 



T V V 

La 0 r 



i=l 

n - 1 



where V 0 is the residual at the center, V r is the residual at a point which is at a 
distance r from the center and n is the number of points. There are several 
methods of determining C(r) of which the concentric circle approach is the 
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simplest. The coefficients B and k can then be computed from the computed 
covariance. 

The expected covariance between the signals at pairs of interpolated stations 
C p can be obtained by substituting the distance between the interpolated stations 

and control stations into the covariance function. 

The expected covariance between all points of mixed observed and 
interpolated signals can be obtained from the covariance function for each 

distance from the interpolated stations to the control stations. 

In this thesis two cases were studied in solving for the geoid height model. 

Case 1. C q = I and = 0 

Case 2. C q ^ I and Cpq 0 

I start with assuming 

C = P 1 = I 
q 

C = 0 
pq 

where P is the weight matrix and I is the unit matrix. 

Then 

X 0 =(A T P A)' 1 A T P x 
and 

S = 0*P (x - A X ) = 0 

p x (r 

Using X to compute residuals, v, and then estimating C^, C pq and P, it is then 
found that 

X, =(A t C' 1 A)‘ 1 A t C‘ 1 x 
i q q 

and 

S = C C 4 (x - A X,) 
p pq q t 
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IV. DATA COLLECTION AND INSTRUMENTATION 



A. PLANNING 

Nine temporary bench marks were established on the NPS campus (Figure 4). 
Bench marks GH7 and GH8, designated as check marks, were established in the 
center of the levelling loop. H of bench marks was obtained by differential 
levelling. Ah was measured with a pair of GPS receivers. 




Figure 4. The temporary bench marks on the NPS campus. 
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Eight permanent bench marks were recovered in the study area (Figure 5) 
[Vertical Control Data, 1961]. Bench mark GWM 27, designated as a check mark, 
was roughly in the center of the survey area for the permanent bench marks. H was 
obtained from the published elevations. Ah was also measured with a pair of GPS 
receivers. 
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B. DATA COLLECTION 
1. Obtaining H 

a. Precise Levelling 

Precise levelling was run on the NPS campus on February 19, 1988, 
and March 24, 1988. The elevation of station TREE was assumed to be zero 
above the MSL, so the orthometric height differences, AH, for each temporary 
bench mark could be observed. The differences of the relative elevations, AH, 
for the forward and the backward sights on the NPS campus are listed in Table 1 . 



Table 1. AH FOR THE NPS CAMPUS 



Bench 


mark 


AH (m) 


From 


To 


Forward 


Backward 


Mean 


TREE 


GH1 


-0.808 


0.808 


-0.808 


GH1 


GH2 


5.388 


-5.384 


5.386 


GH2 


GH3 


3.150 


-3.149 


3.149 


GH3 


GH4 


1.535 


-1.536 


1.535 


GH4 


GH5 


-1.016 


1.016 


-1.016 


GH5 


GH6 


-4.214 


4.215 


-4.214 


GH6 


GH7 


0.489 


-0.489 


0.489 


GH7 


TREE 


-4.031 


4.033 


-4.032 


GH6 


GH2 


0.547 


-0.546 


0.547 


GH7 


GH2 


0.057 


-0.057 


0.057 


GH3 


GH8 


0.716 


-0.718 


0.717 


GH8 


GH4 


0.817 


0.817 


0.817 



To avoid an obstruction near bench mark D 697, which is 10 ft north 
of a 12-inch cedar tree, the temporary bench mark TEMP 1 was established 
nearby. Offset levelling was run on April 8, 1988. The elevation of TEMP 1 
was offset by -0.401 m from the elevation of D 697. 

F 813 which was set in the west end of the south abutment of a steel 
bridge over the Salinas River was offset to temporary bench mark TEMP 2. The 
offset levelling was run on April 9, 1988. The elevation of TEMP 2 was offset 
by -2.218 m from the elevation of F 813. 
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b. Published Elevation 

H values for the larger of our two study areas in and near Monterey, 
California, areas were obtained from the published elevations printed by U.S. 
Department of Commerce Coast and Geodetic Survey Washington D.C. [1961]. 
The geodetic datum used in this publication was the NGVD 29. First-order spirit 
levelling has extended this datum over most of the continent. Although first- 
order lines may be 300 km apart in some western areas, most points in the 
country are no more than 50 km from an estimated first-order bench mark. A 
readjustment of the this network is the NAVD 88. This new adjustment will be 
made to the geopotential surface rather than to sea level, and it will place the 
existing vertical data in a form that makes it most consistent and accessible to the 
user. Changes to older published elevations are not expected to exceed 15 
decimeters [NASA, 1978]. Table 2 lists the orthometric heights of the permanent 
bench marks we used. 



Table 2. H FOR PERMANENT BENCH MARKS UESD IN THIS 
STUDY 



Bench mark 


H(m) 


K 152 


16.965 


B 21 


5.787 


S 812 


15.978 


P 812 


13.134 


D 697 


30.057 


J 697 


85.061 


F 813 


7.983 


GWM27 


63.487 



2. Obtaining Ah 

a. Satellite Observation Plan 

Due to the positions of the satellite orbits during this study, the 
observing window was between 2100 and 0200 hours Pacific Standard Time 
(PST=UTC+8 hours.). The time period was between February 2, 1988, and May 
12, 1988. 
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b. Satellite selection 

The same five satellites (SV) (6, 9, 11, 12 and 13) were used for all 
observations. 

c. Position Dilution of Precision (PDOP) 

The accuracy to which positions are determined using GPS depends 
on two factors: (i) satellite configuration geometry, and (ii) measurement 
accuracy. GPS measurement accuracy represents the combined effect of 
ephemeris uncertainties, propagation errors, clock and timing errors, and 
receiver noise. 

The effect of satellite configuration geometry is expressed by the 
dilution of precision (DOP) factor, which is the ratio of the positioning accuracy 
to the measurement accuracy [Wells, 1987]. 

g = DOP • o 0 



where 

Gq : is the measurement accuracy (standard deviation) 

G : is the positioning accuracy (standard deviation in one coordinate) 
The value of GDOP itself is a composite measure that reflects the 
influence of satellite geometry on the combined accuracy of the estimation of 
observation time (user clock offset) and receiver position [Milliken, 1986]. 

GDOP 2 = PDOP 2 + TDOP 2 

TDOP is the Time Dilution Of Precision, the error in the clock of the receiver 
bias multiplied by the velocity of light. The four best satellites selected by the 
receiver are those with the lowest GDOP. Trimble recommends that the rapidly 
changing PDOP provides better geometry for phase differencing techniques. 
Low or constant PDOP provides weaker solutions. The 4000SX receiver does 
not record GDOP, but it does record PDOP every five minutes. PDOP was 
about 4.5 for all GPS observations. 
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C. INSTRUMENTATION 

1. Levelling 

A Zeiss model Ni-2 level (Ser. # 82377) was used at the temporary 
bench marks for the third-order, class I levelling. It has a 32-diameter 
magnification, produces an erect image, and has stadia constants of 333 or 100. 
A Peg test was performed before levelling. The level error, or c-value was also 
checked before the beginning the levelling. The c-value was -0.006 mm/m which 
was less than +0.05 mm/m, so it was not necessary to adjust the level [Bodnar, 
1975]. The level contains a bubble tube to permit positioning parallel to the 
geoid. When properly set up at a point, the telescope is locked so that it will 
rotate through a 360° arc in a horizontal plane. With the level locked in position 
readings are made on two calibrated staffs held in upright positions ahead of and 
behind the instrument. The difference between readings is the difference in 
elevation between points. Dietzgen Model # 6450 metric rods were used in the 
levelling. These rods are graduated in centimeters. The actual reading is 
estimated to the nearest 1 mm. Rod levels were used to indicate when were 
vertical. 

2. GPS 

a. Trimble 4000SX Receiver 

A complete description of the Trimble 4000SX receiver is given by 
Trimble Navigation [Trimble, 1987a]. NPS operates three Trimble 4000SX GPS 
receivers. For this study I fixed one antenna to the roof of Building 224 on the 
NPS campus, and one was carried to the field. The 4000SX is capable of 
observing the C/A code, integrated Doppler, and carrier beat phases of up to five 
satellites simultaneously. Its ability to use the C/A code allows the receiver to be 
used as a stand alone navigation system to determine positions using Doppler- 
smoothed pseudoranges and velocities [Ashjaee, 1985]. The receiver uses the 
C/A code in a time transfer mode to determine the offset and drift of its own 
clock and thus provide accurate time tags for the observations without the 
requirement of an external atomic clock or synchronization with the receiver at 
the other end of the baseline. The reference position (the geodetic coordinates of 
the antenna) and the practical options chosen must be entered into the receiver via 
the receiver key pad. The 4000SX requires 115 V AC power or 12 V DC 
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power supply. 115 V AC power was used for the NPS campus measurements. 
12 V DC, supplied by a car battery, was used in the field. 

b. Grid Personal Computer (PC) 

For precise relative positioning the 4000SX receiver transmits data 
through an RS-232 port to a microcomputer (Grid PC) for storage on floppy 
disks for post-processing. The Grid PC uses either 115 V AC or 12 V DC. On 
the NPS campus a 115 V AC power supply was used, so measurements could be 
made for four hours. The 12-volt battery in the Grid PC lasts about 120 
minutes, so measurements were taken for only 100 minutes in the field. 

c. Antennas 

Multipath-resistant Trimble microstrip antennas were installed on 
Building 224 on the NPS campus and over the various bench marks. The antenna 
heights from the center of bench marks to the edge of the antenna's ground plane 
were measured before and after GPS observations. The field antenna was 
mounted on a tripod with a tribach and optical plummet for centering and 
levelling the antenna. Arrows on the antennas were oriented to the north at both 
stations using a magnetic compass. 

d. Meteorological Instruments 

A barometer and a sling psychrometer were used to measure 
atmospheric pressure, relative humidity and air temperature at each field site. 

e. The steel tape 

A three-meter steel tape was used to measure the antenna height. 
This tape can be read to 1 cm. Readings are estimated to 1 mm. 

D. SOFTWARE 

A complete description of the Trimble-supplied Trimvec software is to be 
found in Trimble Navigation [Trimble, 1987b]. The software provides data 
logging, baseline computation and datum transformation programs. These 
operate with an IBM compatible personal computers. A printer and a hard disk 
are used with the planning and processing programs. The following programs 
are used on 3-1/2 inch micro-floppies: 

1. Data Logger on Disk 1 

In relative positioning, the receiver is controlled from the Grid PC by 
version D of Trimble’s Gridlog5.bat program. Each observation session was 
initialized to log data when a minimum of four satellites were 15° above the 
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antenna’s horizon. Five satellites were designated for each observing session. 
The observables and receiver clock parameters were logged to a floppy disk 
every 15 seconds and the C/A code-determined antenna position and PDOP every 
five minutes. The GPS navigation message was logged to a separate file at the 
beginning of the session. 

2. Post-processor on Disk 2 

The data was processed using the Trimvec Trim640 program, Revision 
AB. Trim640 is a relative positioning, post-processing program that provides 
triple and double difference solutions for two sets of the 4000SX carrier phase 
data logged simultaneously at two stations. Trim640 adopts the best C/A code 
position during the data loading. Only the broadcast ephemeris can be used to 
compute fixed orbit satellite positions. Trim640 limits processing to 700 
epoches, so the first 700 epoches for each observing session are used. The 
antenna on the roof of Building 224 on the NPS campus was used as reference 
station and its coordinates were kept fixed. Permanent bench marks were chosen 
as Trim640 reference stations for the off campus area because the observing 
sessions at them were shorter than the observing session at Building 224 on the 
NPS campus. 

3. Geoid.exe on Disk 3 

This program computes global geoid height, N 0 (X,Y,Z), at any point 

with an accuracy of a few meters [Trimble, 1987b]. The global geoid heights, 
N 0 (X,Y,Z), are based on Rapp's 1978 360 x 360 model, an harmonic expansion 

of the geopotential referred to GRS 80. Heights are suitable for showing relative 
shape and trends in geoid. Global geoid height for the Monterey Bay area from 
36° 34’ N to 36° 49’ N latitude, and from 121° 34' W to 121° 55' W longitude is 
given in Figure 6. A Fortran program GPSCON (Appendix A) was used for the 
contoured plot. 

A second program also named Geoid.exe was obtained from the NGS of 
NOAA's National Ocean Service Charting and Geodetic Services uses the same 
Rapp’s 1986 360 x 360 model. The global geoid height in the Monterey Bay area 
from the NGS program is shown in Figure 7. Because the NGS program rounds 
off to the nearest decimeters, there are some steps on the curves. The differences 
between Figure 6 and 7 may be attributed to differences between Rapp’s 1978 
and 1986 models. 
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Figure 6. The global geoid height in the Monterey Bay area calculated 
using the Trimvec Geoid.exe program. 




Figure 7. The global geoid height in the Monterey Bay area calculated 
using the NGS Geoid.exe program. 
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4. Satellite Visibility Program on Disk 3 

The program Stvis.exe is a system for generating visibility data for the 
GPS satellites. The program takes as input the observer’s latitude, longitude, a 
mnemonic name for the location, data for calculations and optionally the time zone 
offset from Universal Time Coordinates (UTC). NPS in Monterey, California, was 
used as the reference station ( 36° 35' 56.1" N, 121° 52' 36.0” W, and altitude -19 
m). Satellite orbit data is contained in the file, Almanac.dat, which is produced by 
processing data collected from the actual satellites. The satellite visibility chart 
shows the tracks of the GPS satellites during the planned observation session 
(Figure 8). 




Figure 8. Sky plots of satellite tracks for the Monterey Bay area. 
Elevation angles are dotted concentric circles. Zenith is at the center. 
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V. DATA PROCESSING 



A. ORTHOMETRIC HEIGHT ON NPS CAMPUS 

In calculating the most probable elevations for the temporary bench marks, 
station TREE was used as a reference, where the elevation was assumed to be 0 
meters above MSL. The interlocking levelling circuit on the NPS campus is 
shown in Figure 9. 




Figure 9. The interlocking levelling circuit on the NPS campus. 



A constrained least squares adjustment was used to do the computation. A 
weight of 100 was assumed at the TREE and weights of 1 were used for all 
observations. The 13 observation equations are: 
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The observation equations in matrix form are [Wolf, 1980]: 

P A H = PL + PV 

The weight matrix P 13xl3 , coefficient matrix A 13x9 , observation matrix L 13xl , 
and residual matrix V 13xl are: 
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v J 



L = 



/ 0.000 \ 






-0.808 




V„ 


5.386 




v 3 

\T 


3.149 




V 4 


1.535 




V 5 


-1.016 




V 6 


-4.214 


V = 


V 7 


0.489 




V 8 


-4.032 




v 0 


0.547 




9 

V „ 


0.057 




V 10 
\T 


0.717 




V 11 


V0. 817 




^12 



The program Lobs.basic [Jeyapalan, 1988] was used to find H for the temporary 
bench marks. The program (Appendix B) gives the H and V matrices: 



32 







/ 0.00000 \ 
0.00007 


/ 0.000 \ 




0.00007 


-0.808 1 




-0.00022 


4.579 




-0.00044 


7.728 




-0.00022 


9.262 


V = 


-0.00022 


8.246 




0.00007 


4.032 




-0.00019 


4.521 




-0.00010 


V 8.445 J 




-0.00010 
0.00022 
\0. 00023 J 



The standard deviation of the observation equations is 0.36 mm. 

B. GPS DATA PROCESSING 

1. Automatic Processing Mode 

Trimble recommends that an automatic processing mode be used if more 
than one hour of data from four or five GPS satellites is taken at each site for 
baselines lengths up to 30 km. 

The output files contain triple difference, double difference float and 
double difference fixed solutions. The double difference fixed solution provides 
the most accurate solution if the following conditions are met: 

a. Integer bias search indicates the ratio sum-of-squares must be greater than 
3.0; 

b. RMS (cycles) less than 0.05; and 

c. Difference between the float and fixed solutions is less than 10 cm, in any 
component of the baseline (X,Y,or Z). 

In the GPS data processing all the solutions from the bench marks 
satisfied these conditions except station J 697 for which RMS was 0.068. The 
distance from NPS Building 224 to J 697 is 33 km. Trimble recommends a 
baseline length greater than 30 km when the fixed solution does not meet the 
above conditions. The float solution should be used provided the RMS of fit is 
better than 0.08. Ah of double difference fixed solutions were used for the bench 
marks in this study. The double difference fixed solutions are shown in Table 3. 
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Table 3. DOUBLE DIFFERENCE FIXED SOLUTIONS 



Bench 

mark 


Slope distance 
from Bldg. 224 
(m)__ 


RMS 

(cycles) 


Ratio 


Coordinates difference between 
fixed and float solution 
AX (m) AY (m) AZ (m) 


TREE 


66 


0.039 


20.6 


-2.0 


0.3 


-1.4 


GH1 


98 


0.014 


257.3 


-0.5 


0.4 


0.0 


GH2 


256 


0.025 


94.7 


-3.2 


3.7 


0.3 


GH3 


536 


0.014 


276.4 


-1.1 


0.6 


-0.1 


GH4 


547 


0.019 


151.0 


-0.7 


0.4 


-0.2 


GH5 


424 


0.014 


245.7 


0.0 


-1.6 


-1.2 


GH6 


166 


0.015 


254.9 


-0.7 


0.5 


-0.1 


GH7 


169 


0.015 


288.6 


-0.7 


0.3 


-0.1 


GH8 


413 


0.016 


160.0 


-0.4 


0.5 


0.1 


K 152 


14914 


0.035 


9.7 


-0.3 


0.1 


-0.4 


B 21 


2082 


0.014 


193.5 


0.4 


-0.2 


0.1 


S 812 


17062 


0.020 


15.8 


-1.2 


1.7 


0.1 


P 812 


20011 


0.020 


31.7 


0.9 


-1.4 


0.4 


D 697 


25249 


0.045 


11.1 


8.2 


9.3 


0.4 


J 697 


32692 


0.068 


6.3 


13.8 


-10.7 


2.6 


F 813 


16861 


0.022 


35.9 


7.1 


-5.2 


0.1 


GWM27 


16687 


0.049 


15.3 


4.4 


-5.7 


0.5 



Default parameters are assumed for automatic processing. These 
parameters have a minimum elevation mask of 15° for double differences. 
Station one coordinates were the best code positions during data logging. Actual 
surface meteorological values were used in the processing. 

For sets of data covering more than one hour, epoch increments of five 
or ten provide full precision. Automatic processing begins with several iterations 
of triple differences to establish a starting value for the double difference 
solution. A cycle slip fixer processes every epoch. 

Double difference solutions begin with several iterations of the float 
solution where biases, as well as baseline components, are computed. It is called 
the float solution, since the biases are allowed to float and are not constrained to 
be integers. Each iteration should indicate convergence toward zero for 
observed minus computed phases. The default elevation mask is 15° for 
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processing double differences. After completion of the float solution, 
correlations and biases are computed. The integer search algorithm begins to 
determine the proper integer values for the bias. After the biases are set to their 
integer values, the processing is repeated to compute the fixed solution. Ellipsoid 
height differences Ah for the temporary bench marks with respect to NPS 
Building 224, found using double difference fixed solutions is given in Table 4. 



Table 4. Ah AND WGS 84 COORDINATES FOR THE 
TEMPORARY BENCH MARKS. 



Bench mark 


Ah (m) 


X (m) 


Y(m) 


Z (m) 


TREE 


-8.3532 


-2707298.015 


-4353450.286 


3781781.432 


GH1 


-9.1814 


-2707385.942 


-4353407.659 


3781790.139 


GH2 


- 3.7766 


-2707502.139 


-4353548.202 


3781549.491 


GH3 


-0.6382 


-2707528.391 


-4353750.407 


3781314.047 


GH4 


0.9055 


-2707367.536 


-4353813.483 


3781306.605 


GH5 


-0.1099 


-2707202.581 


-4353799.531 


3781493.061 


GH6 


-4.3397 


-2707252.929 


-4353587.524 


3781666.956 


GH7 


-3.8509 


-2707406.088 


-4353553.901 


3781589.428 


GH8 


0.0843 


-2707329.452 


-4353759.387 


3781428.906 



Ah for the permanent bench marks with respect to NPS Building 224 
using double difference fixed solutions are given in Table 5. 
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Table 5. Ah AND WGS 84 COORDINATES FOR THE 
PERMANENT BENCH MARKS 



Bench mark 


Ah (m) 


X (m) 


Y(m) 


Z(m) 


K 152 


2.2398 


-2696890.306 


-4350931.132 


3792065.393 


B21 


-8.4275 


-2708436.436 


-4351973.362 


3782674.213 


S 812 


1.6335 


-2692137.111 


-4360870.528 


3784094.277 


P 812 


-1.3569 


-2689245.909 


-4360679.825 


3786312.651 


D 697 


15.1728 


-2685153.665 


-4357222.820 


3793175.171 


J 697 


71.3023 


-2679276.491 


-4356508.580 


3798216.777 


F 813 


-6.7070 


-2695613.779 


-4350454.843 


3793463.222 


GWM 27 


49.1051 


-2692408.721 


-4360427.654 


3784433.846 



2. Obtaining Geoid Height 

Geoid heights of bench marks were obtained by using the Geoid.exe 
program from NGS or Trimvec (Table 6). 



Table 6. GEOID HEIGHTS 



Bench mark 


NGS (m) 


Trimvec (m) 


TREE 


-34.7 


-34.30428 


GH1 


-34.7 


-34.30671 


GH2 


-34.7 


-34.31508 


GH3 


-34.7 


-34.32072 


GH4 


- 34.7 


-34.31605 


GH5 


-34.7 


-34.30745 


GH6 


-34.7 


-34.30532 


GH7 


-34.7 


-34.31140 


GH8 


-34.7 


-34.31247 


K 152 


-34.1 


-33.83303 


B 21 


-34.7 


-34.31918 


S 812 


-33.9 


- 33.86473 


P 812 


-33.8 


- 33.76630 


D 697 


-33.6 


-33.58265 


J 697 


-33.3 


-33.42047 


F 813 


-34.0 


- 33.77899 


GWM 27 


-33.9 


-33.86547 
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3. Determination of Local Geoid Model 

The method of collocation was used to solve the local geoid model. The 
procedure is described as follows : 

a. Reference Station Selection 

Station TREE was selected as a reference station on the NPS campus. 
Coordinate differences AX, AY and AZ in WGS 84, for the temporary bench 
mark with respect to the TREE were computed. These are shown in Table 7. 



Table 7. AX, AY AND AZ ON THE NPS CAMPUS 



Bench mark 


AX (m) 


AY (m) 


AZ (m) 


TREE 


0.000 


0.000 


0.000 


GH1 


-87.927 


42.627 


8.707 


GH2 


-204.124 


-97.916 


-231.941 


GH3 


-230.376 


-300.121 


-467.385 


GH4 


-69.521 


-363.197 


-474.827 


GH5 


95.434 


-349.245 


-288.371 


GH6 


45.085 


-137.238 


-114.476 


GH7 


-108.073 


-103.615 


-192.004 


GH8 


-31.437 


-309.101 


-352.526 



Station K 152 was selected as a reference station in the large off 
campus area. The coordinate differences, AX, AY and AZ between the 
permanent bench marks and K 152 were computed in WGS 84 (Table 8). 
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Table 8. AX, AY AND AZ IN THE OFF CAMPUS AREA 



Bench mark 


AX (m) 


AY (m) 


AZ (m) 


K 152 


0.000 


0.000 


0.000 


B 21 


-11546.130 


-1042.230 


-9391.180 


S 812 


4753.195 


-9939.396 


-7971.116 


P 812 


7644.397 


- 9748.693 


-5752.742 


D 697 


11736.641 


-6291.688 


1109.778 


J 697 


17613.815 


-5577.448 


6151.384 


F 813 


1276.527 


476.289 


1397.829 


GWM27 


4481.585 


-9496.522 


-7631.547 



b. Determine the Local Geoid Model 

I start with assuming the signal, S, equals zero and the weight matrix 
P, equals the unit matrix. Seven control marks in both areas led to seven 
observation equations. The coordinate differences between the reference marks, 
which are TREE on the NPS campus and K 152 in the off campus area, are the 
coefficients of the parameters. The observation equation is given by 

AN = H + N q (X,Y,Z) - Ah 



where 

H : the orthometric height, from levelling 

N 0 (X,Y,Z) : the global geoid height, obtained form NGS or Trimvec 
Ah : the ellipsoid height difference with respect to the reference station 
The local variation of the geoid height, AN, is also given by 

AN = h Q + combination of coordinate differences 

h 0 , the ellipsoid height, includes a constant correction to N 0 (X,Y,Z) at the 
reference station. AN for temporary bench marks, N 0 (X,Y,Z) from the NGS 
(designated as NAN) and N 0 (X,Y,Z) from the Trimvec (designated as TAN) are 
given in Table 9. 
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Table 9. AN ON THE NPS CAMPUS 



Bench mark 


NAN (m) 


TAN (m) 


TREE 


-26.3468 


-25.9511 


GH1 


-26.3266 


-25.9333 


GH2 


-26.3444 


-25.9595 


GH3 


-26.3338 


-25.9545 


GH4 


-26.3435 


-25.9595 


GH5 


-26.3441 


-25.9515 


GH6 


-26.3283 


-25.9336 



AN for permanent bench marks with N 0 (X,Y,Z) from the NGS (designated as 
NAN) and N 0 (X,Y,Z) from the Trimvec (designated as TAN) are given in Table 
10 . 



Table 10. AN IN THE OFF CAMPUS AREA 



Bench mark 


NAN (m) 


TAN (m) 


K 152 


-19.3748 


-19.1078 


B 21 


-20.4855 


-20.1047 


S 812 


-19.5555 


-19.5202 


P 812 


-19.3091 


-19.2754 


D 697 


-18.7158 


-18.6984 


J 697 


-19.5413 


-19.6618 


F 813 


-19.3100 


-19.0890 



A least squares adjustment program Lobs.basic was used to determine 
the geoid model which has the smallest standard deviation. This was done by 
using the different combinations of the coordinate differences. GPS data from 
the off campus area were studied to determine the local geoid model. N 0 (X,Y,Z) 

from the Trimvec program were used to compute AN. 
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Four-parameter combinations and their standard deviations for the off 
campus area are shown in Table 11. 

Table 11. FOUR PARAMETERS 



Parameters 


Standard deviation (m) 


h Q + aAX + bAY + cAZ 


0.45 


h 0 + aAX + bAY + cAY 2 


0.52 


h 0 + aAX + bAY + cAX 2 


0.34 


h 0 + aAX + bAY + cAXAY 


0.47 


h Q + aAX + bAZ + cAY 2 


0.52 


h Q + aAX + bAZ + cAX 2 


0.34 


h Q + aAX + bAZ + cAXAY 


0.47 


h 0 + aAX + bAY 2 + cAX 2 


0.31 


h 0 + aAX + bAX 2 + cAXAY 


0.37 


h 0 + aAY + bAZ + cAY 2 


0.52 


h Q + aAY + bAZ + cAX 2 


0.34 


h Q + aAY + bAZ + cAXAY 


0.46 


h 0 + aAY + bAY 2 + cAX 2 


0.43 


h Q + aAY+ bAY 2 + cAXAY 


0.58 


h Q + aAY+ bAX 2 + cAXAY 


0.28 


h Q + aAZ + bAY 2 + cAX 2 


0.37 


h 0 + aAZ + bAX 2 + cAXAY 


0.32 


h Q + aAZ + bAY 2 + cAXAY 


0.49 


h Q + aAY 2 + bAX 2 + cAXAY 


0.20 * 



* The smallest standard deviation. 
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Five-parameter combinations and their standard deviations for the off campus 
area are given in Table 12. 



Table 12. FIVE PARAMETERS 



Parameters 


Standard deviation (m) 


h Q + aAX + bAY + cAZ + dAX 2 


0.43 


h Q + aAX + bAY + cAZ + dAY 2 


0.57 


h Q + aAX + bAY + cAZ + dAXAY 


0.39 


h Q + aAX + bAZ + cAX 2 + dAY 2 


0.57 


h Q + aAX + bAZ + cAY 2 + dAXAY 


0.59 


h 0 + aAX + bAX 2 + cAY 2 + dAXAY 


0.22 


h Q + aAX + bAZ + cAX 2 + dAY 2 


0.57 


h Q + aAY + bAZ + cAX 2 + dAXAY 


0.34 


h Q + aAY + bAZ + cAX 2 + dAXAY 


0.34 


h Q + aAY + bAZ + cAX 2 + dAY 2 


0.56 


h Q + aAY + bAX 2 + cAY 2 + dAXAY 


0.10 * 


h Q + aAZ + bAX 2 + cAY 2 + dAXAY 


0.19 



* The smallest standard deviation. 



Six-parameter combinations and their standard deviations for the off campus area 
are given in Table 13. 



Table 13. SIX PARAMETERS 



Parameters 


Standard deviation (m) 


h Q + aAX + bAY + cAZ + dAX 2 + eAY 2 


0.26 


h Q + aAX + bAY + cAZ + dAX 2 + eAXAY 


0.59 


h Q + aAX + bAY + cAX 2 + dAY 2 + eAXAY 


0.13 


h Q + aAX + bAZ + cAX 2 + dAY 2 + eAXAY 


0.12 * 


h Q + aAY + bAZ + cAX 2 + dAY 2 + eAXAY 


0.14 



* The smallest standard deviation. 
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For seven observation equations and seven parameters the degree of freedom 
equals zero, so the standard deviation is undefined. Seven-parameter 
combinations and their standard deviations for the off campus area are given in 
Table 14. 



Table 14. SEVEN PARAMETERS 



Parameters 


Standard deviation (m) 


h Q + aAX + bAY + cAZ + dAX 2 + eAY 2 + fAXAY 


undefined 



The best standard error of the four-parameter geoid model is greater 
than for five-parameter or six-parameter geoid models. The seven-parameter 
geoid model did not converge, so five-parameter and six-parameter geoid models 
were chosen for this study. The five-parameter geoid model is given by 

N = N q (X,Y,Z) + h 0 + aAY + bAX 2 + cAY 2 + dAXAY 

and the six-parameter geoid height is given by 

N = N 0 (X,Y,Z) + h Q + a’AX + b’AZ + c'AX 2 + d'AY 2 + e’AXAY 

c. Evaluating the Geoid Model Using Check Marks 

In order to evaluate the accuracy of the geoid model, the geoid height 
model can be rearranged to calculate the orthometric height of check marks 
which are GH7 and GH8 on the NPS campus and GWM 27 in the off campus 
area. 



H = h 0 + Ah - N 0 (X,Y,Z) + aAY + bAX 2 + cAY 2 + dAXAY 
or 

H = h 0 + Ah - N 0 (X,Y,Z) + a’AX + b'AZ + c’AX 2 + d’AY 2 + e’AXAY 

The h 0 , a, b, c and d of the five-parameter geoid model, and h 0 , a', b', c', d' and 

e' of the six-parameter geoid model, can be solved by least squares adjustments. 
The results for h 0 , a, b, c and d for the five-parameter model using the seven 
control marks in the off campus area are given in Table 15. The results for h 0 , 
a, b, c and d are shown in the NGS column, where the N 0 (X,Y,Z) was from the 
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NGS Geoid.exe program. The results for h 0 , a, b, c, d and standard deviation, a, 
are shown in the Trimvec column, where the N 0 (X,Y,Z) was from the Trimvec 
Geoid.exe program. 



Table 15. h 0 , a, b, c, d, a FOR SEVEN CONTROL MARKS 



Parameter 


NGS (m) 


Trimvec (m) 


ho 


-19.26895 


-19.00915 


a 


-2.344959E-04 


-2.736598E-04 


b 


-8.523017E-09 


-8.528133E-09 


c 


-3.5743 14E-08 


-3.905052E-08 


d 


-2.211550E-08 


-1.728313E-08 


a 


0.130518 


0.095723 



The results of h 0 , a, b, c, d and a of the five-parameter geoid model 
using the six control marks (excluding B 21) in the off campus area are given in 
Table 16. 



Table 16. h 0 , a, b, c, d, a FOR SIX CONTROL MARKS 



Parameter 


NGS (m) 


Trimvec (m) 


ho 


-19.25415 


-19.01411 


a 


-2.752543E-04 


-2.600850E-04 


b 


-7.735253E-09 


-8.792540E-09 


c 


-3.777950E-08 


-3.836067E-08 


d 


-1.798617E-08 


-1.867193E-08 


CT 


0.171919 


0.133516 



The results of h 0 , a, b, c, d and a of the five-parameter geoid model 
using the six control marks (excluding J 697) in the off campus area are given in 
Table 17. 
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Table 17. h 0 , a, b, c, d, a FOR SIX CONTROL MARKS 



Parameter 


NGS (m) 


Trimvec (m) 


h 0 


-19.265355 


-19.03487 


a 


-2.503395E-04 


-1.977533E-04 


b 


-8.731945E-09 


-7.52641 0E-09 


c 


-3.711466E-08 


-3.246532E-08 


d 


-2.179058E-08 


-1. 88279 IE-08 


a 


0.1841317 


0.1206984 



The results of h 0 , a, b, c and d of the five-parameter geoid model using 
the five control marks (excluding B 21 and J 697) in the off campus area are 
given in Table 18. The standard deviation is undefined since degree of freedom 
equals zero. 



Table 18. h 0 , a, b, c, d FOR FIVE CONTROL MARKS 



Parameter 


NGS (m) 


Trimvec (m) 


h 0 


-19.37497 


-19.10791 


a 


1.032352E-04 


3.397465E-05 


b 


8.1 19969E-09 


3.521564E-09 


c 


7.377821E-09 


-3.339665E-09 


d 


1.338776E-09 


-3.667083E-09 



The results of h 0 , a, b, c, d and a of the five-parameter geoid model 
using the seven control marks on the NPS campus are given in Table 19. 



Table 19. h 0 , a, b, c, d, a FOR SEVEN CONTROL MARKS 



Parameter 


NGS (m) 


Trimvec (m) 


h 0 


-26.33543 


-25.94105 


a 


2.976507E-06 


1.192093E-07 


b 


-5.419133E-08 


-1.763692E-07 


c 


-3.601599E-08 


-8.489588E-08 


d 


6.30971 IE-08 


-3.702007E-08 


c 


0.013773 


0.014612 
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The results of h 0 , a', b', c', d', e* and o of the six-parameter geoid model 
using the seven control marks in the off campus area are shown in Table 20. 



Table 20. h 0 , a', b\ c’, d', e\ a FOR SEVEN CONTROL MARKS 



Parameter 


NGS (m) 


Trimvec (m) 


h o 


- 19.27388 


-19.01660 


a’ 


2.083 182E-04 


1.842380E-04 


b’ 


-2.5615 10E-04 


-2.520234E-04 


c' 


- 7.731387E-09 


-8.335974E-09 


cl 1 


-3.767036E-08 


-3.960304E-08 


e’ 


-1.249646E-08 


-1.521767E-08 


a 


0.137149 


0.123960 



The results of h 0 , a', b\ c\ d' and e' of the six-parameter geoid model 
using the six control marks (excluding B 21) in the off campus area are given in 
Table 21. The standard deviation is undefined since degree of freedom equals 
zero. 



Table 21. h 0 , a*, b’, c’, d’, e' FOR SIX CONTROL MARKS 



Parameter 


NGS (m) 


Trimvec (m) 


h o 


-19.37489 


-19.10785 


a’ 


2.57641 IE-04 


2.288222E-04 


b' 


-1.691282E-04 


-1. 73475 6E-04 


c' 


-1. 08307 2E-08 


-1.1 13654E-08 


d 1 


-2.817797E-08 


-3.1 0274 IE-08 


e' 


-5.824404E-09 


-9.185896E-09 



The results of h 0 , a', b\ c', d' and e' of the six-parameter geoid model 
using the six control marks (excluding J 697) in the off campus area are given in 
Table 22. The standard deviation is undefined since degree of freedom equals 
zero. 
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Table 22. h 0 , a’, b’, c’, d*, e* FOR SIX CONTROL MARKS 



Parameter 


NGS (m) 


Trimvec (m) 


ho 


-19.37466 


-19.10774 


a’ 


1.065433E-04 


9.238720E-05 


b' 


-4.604459E-05 


-6.231666E-05 


c' 


-2.012712E-09 


-3.170499E-09 


d’ 


-1.141234E-08 


-1.586523E-08 


e’ 


-2.528395E-09 


- 6.206392E-09 



The results of h 0 , a', b', c', d', e' and a of the six-parameter geoid model 
using the seven control marks on the NPS campus are given in Table 23. 



Table 23. h 0 , a*, b’, c’, d', e’, a FOR SEVEN CONTROL MARKS 



Parameter 


NGS (m) 


Trimvec (m) 


ho 


-26.33389 


-25.93820 


a’ 


1.645088E-05 


4.556775E-05 


b' 


4.225970E-05 


6.1 6908 IE-05 


c’ 


6.856863E-08 


6.391201E-08 


d’ 


6.00703 IE-08 


5.878974E-08 


e’ 


1.781 154E-07 


1.767185E-07 


a 


0.018859 


0.018872 



Fortran program GPSDIS (Appendix C) was used to calculate the 
orthometric heights of the check marks. The results are discussed in the Chapter 
VII. 

d. Obtaining the Weight Matrix 

Since the standard error of the geoid model is about + 2 cm, which is 
the same or better than GPS Ah accuracy, it is not actually necessary to do 
further computation. Nevertheless, the weight matrices, P, were computed. 

The five-parameter geoid model used for the NPS campus and the 
six-parameter geoid model used for the off campus area were employed in 
obtaining P. The covariance function, C(r), of the control marks (p. 18) can be 
obtained by using the least squares residuals. The coefficients B and k of C(r) 
can be obtained by solving simultaneously the equations for C(r) with different 
distances from the centers. 
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(1) P for NPS Campus . The residuals of the five-parameter geoid 
model of the control marks on the NPS campus are given in Table 24. The 
global geoid heights were from the NGS. 



Table 24. RESIDUALS FOR FIVE-PARAMETER MODEL 



Control mark 


# of residual 


Residual (m) 


TREE 


v, 


1.137352E-02 


GH1 


V 2 


-9.420395E-03 


GH2 


V 3 


7.339478E-03 


GH3 


V 4 


-4.278183E-03 


GH4 


v 5 


3.572464E-03 


GH5 


V 6 


6.446839E-04 


GH6 


V 7 


- 8.714676E-03 



For r = 593 m 



C(r ) = 



V V + V V + V V + V V, + V V, 

14 15 24 25 46 



-2.194272E-03 



For r 2 = 3 1 9 m 



V, V ,+ V , V ,W V+V V+V V +V V +V V.+V V +V 
C(r )= — ^ — 2 — ^ — 2-6 — 2 _j — 4_6 — 5 — ~= 3.94581 3E-6 

2 o 



Thus the equations of covariance can be written as 



CO^) = 0.0138 + B sinCkTj) 



and 



C(r 2 ) = 0.0138 + B sin(kr 2 ) 

Approximate solutions for coefficients B and k can be found by using the 
Maclaurin series expansion. The solved covariance function is given by 

C(r) = 0.0138 - 0.0160 sin(1.7670 r) 
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C q can be obtained by using the C(r) with the distances between the control 

marks. Fortran program DISTCO (Appendix D) was used to compute the 
distances between the control marks, as well as C q and, hence, the weight matrix, 
P = Cq 1 . The Matlab program on the NPS mainframe computer was used to find 
the inverse of C q (P = C q _1 ). P for the five-parameter geoid model for the NPS 
campus is given by 



s' 



P = 



N 



162 



ZEROS 



142 



78 



98 



102 



ZEROS 



87 



100 



\ 



(2) P for off Campus Area . The residuals of the six-parameter geoid 
model of the control marks in the off campus area are given in Table 25. The 
global geoid heights were calculated using the Trimvec program. 



Table 25. RESIDUALS FOR SIX-PARAMETER MODEL 



Control mark 


# of residual 


Residual (m) 


K 152 


v, 


9.120178E-02 


B 21 


V 2 


-9.773254E-03 


S 812 


V 3 


6.391526E-03 


P 812 


V 4 


1.983643E-04 


D 697 


V 5 


-2.77996 IE-02 


J 697 


V 6 


1.685524E-02 


F 813 


V 7 


-7.65171 IE-02 



For Tj = 13994 m 



V V + V , V . + V 1 V . + V , V < + V . V 7 + V 4 V 7 

QTj) = — *—2 L - J 1415 l — L 1- 2 - = -6.659883E-04 
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For r 2 = 19230 m 




vv+vv+vv 

1 6 2 3 3 6 



= 7.912463E-04 



2 



Thus the equations of covariance can be written as 



C(r ) = 0.0124 + B sin(kr ) 



and 



C(r 2 ) = 0.0124 + B sin(kr 2 ) 



Approximate solutions for coefficients B and k can be found by using a 
Maclaurin series expansion. The solved covariance function is given by 

C(r) = 0.0124 - 0.1329 sin(2.8278 r) 

Similarly, C q can be obtained by using the C(r) with the distances between the 

control marks. Fortran program DISTCO (Appendix D) was used to compute 
the distances between the control marks and C q , to obtain the weight matrix (P = 
Cq' 1 ). The Matlab program on the NPS mainframe computer was used to find the 
inverse of C q (P = Cq' 1 ). P for the six-parameter geoid model of the off campus 
area is given by 



25 



ZEROS 



15 



15 



P 



16 



10 



ZEROS 



17 



25 



j 
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VI. EVALUATION OF GPS RECEIVER ERRORS 



Factors affecting the GPS measured ellipsoid height differences, Ah, 
include: (1) comparing GPS differencing solutions, (2) standard error of GPS 
observations, (3) corrections for surface meteorological values, and (4) 
observation durations for GPS. 

A. COMPARING GPS DIFFERENCING SOLUTIONS 

1. Ellipsoid Height Differences Tested at a Fixed Position 

In order to determine the best solution for Ah using GPS, the Ah data 
from a Trimble 4000SX GPS receiver were collected in the K parking lot on the 
NPS campus from 0633 to 1007 UTC on February 19, 1988. The design of the 
experiments is shown in the Figure 10. 




Fixed point on the ground 



Figure 10. Illustration of the Ah test of Trimble 4000SX GPS 
receiver. 
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A plumb bob was used to point to a fixed position (36° 35' 56" N, 121° 
52' 36" W) on the ground. The observations were taken by changing the antenna 
height over a distance of about 60 cm, which was measured with a tape every 
half hour. The ellipsoid height differences, Ah, and the orthometric height 
differences, AH, should agree with the change of antenna height for a fixed 
position. Triple difference solutions, Ah TRI , are given in Table 26. Table also 
gives the differences between the Ah TRI at succeeding observation times, DAh TRI ; 
double difference float solutions, Ah FLT ; the differences between the Ah FLT at 
succeeding observation times, DAh FLT ; and double difference fixed solutions, 
Ah FIX , and the differences between the Ah nx at succeeding observation times, 

DAhpjx- 



Table 26. Ah AND DAh FOR DIFFERENT ANTENNA HEIGHTS 



TIME (UTC) 


Ah-fRi 


DAhyRi 


Ah pur 


DAhpLj 


^FIX 


DAhppx 


(19-2-1988) 


(m) 


(m) 


(m) 


(m) 


(m) 


(m) 


0633 - 0712 


4.4410 


-0.7189 


4.4429 


-0.6339 


4.5045 


-0.5463 


0736 - 0808 


5.1599 


-0.3973 


5.0768 


-0.5676 


5.0508 


-0.5964 


0824 - 0855 


5.5572 


-0.1954 


5.6444 


-0.7144 


5.6472 


-0.7030 


0911 -0943 


5.7526 


-0.4531 


6.3588 


-0.2920 


6.3502 


-0.2660 


0945 - 1007 


6.2057 




6.6508 




6.6162 





The differences between the GPS DAh and the differences of the 
measured ellipsoid height differences, DAh MEA , are given in Table 27. 
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Table 27. THE DIFFERENCE OF Ah 



ANTENNA 
HEIGHT (m) 


DAIimea 

(m) 


DAh MEA - DAh TRI 
(m) 


DAh MEA - DAh^ 
(m) 


DAh MEA - DAhpjx 

(m) 


2.809 










2.262 


-0.5470 


0.1719 


0.0869 


-0.0007 


1.666 


-0.5960 


-0.1987 


-0.0284 


0.0004 


1.074 


-0.5920 


-0.3966 


0.1224 


0.1110 


0.484 


-0.5900 


-0.1369 


-0.2980 


-0.3240 



The results show that the double difference fixed solution is the best. 
The difference between the measured ellipsoid height difference and the GPS 
ellipsoid height difference was about 1 mm. The last two results were not good 
because they had multipath effects and imaging effects from the ground (antenna 
heights 1.074 m and 0.484 m), so the solutions of AH compared to the true 
differences were large. This indicates that the antenna should be at least 1 m 
above the ground. 

2. Ellipsoid Height Differences Tested at Different Positions 

Since the geoid on the NPS campus has a flat geoid slope, the differences 
between Ah, DAh, and the orthometric height difference, AH, of the temporary 
bench marks can be neglected. The Ah and DAh from the GPS observations are 
given in Table 28. The differences between the GPS DAh and the AH are given 
in Table 29. 
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Table 28. Ah AND DAh ON THE NPS CAMPUS 



Bench 


Ah T Ri 


DAhiRi 


AhpLx 


DAhpLT 


Ah fix 


DAhpp,. 


mark 


(m) 


(m) 


(m) 


(m) 


(m) 


(m) 


TREE 


-8.3358 


0.8501 


-8.3375 


0.8313 


-8.3532 


0.8282 


GH1 


-9.1859 


-5.3472 


-9.1688 


-5.3690 


-9.1814 


-5.4048 


GH2 


-3.8387 


-3.1994 


-3.7998 


-3.1624 


-3.7766 


-3.1384 


GH3 


-0.6393 


-1.5241 


-0.6374 


-1.5441 


-0.6382 


-1.5437 


GH4 


0.8848 


0.9777 


0.9067 


1.0171 


0.9055 


1.0154 


GH5 


-0.0929 


4.2655 


-0.1104 


4.2303 


-0.1099 


4.2298 


GH6 


-4.3584 


-0.4932 


-4.3407 


-0.4916 


-4.3397 


-0.4888 


GH7 


-3.8652 


-3.9412 


-3.8491 


-3.9360 


-3.8509 


-3.9352 


GH8 


0.0760 




0.0869 




0.0843 





Table 29. DIFFERENCES BETWEEN DAh AND AH 



Bench 


H 


AH 


AH - DAhyRj 


AH - DAh^ 


AH - DAhp^ 


mark 


(m) 


(m) 


(m) 


(m) 


(m) 


TREE 


0.000 


0.808 


-0.0423 


-0.0235 


-0.0204 


GH1 


-0.808 


-5.387 


-0.0398 


-0.0180 


0.0178 


GH2 


4.579 


-3.149 


0.0504 


0.0134 


-0.0106 


GH3 


7.728 


-1.534 


-0.0099 


0.0101 


0.0097 


GH4 


9.262 


1.016 


0.0383 


-0.0011 


0.0006 


GH5 


8.246 


4.214 


-0.0515 


-0.0163 


-0.0158 


GH6 


4.032 


-0.489 


0.0042 


0.0026 


-0.0002 


GH7 


4.521 


-3.924 


0.0172 


0.0120 


0.0112 


GH8 


8.445 











53 



It is clear that the double difference fixed solutions give the best solution 
in a flat geoid slope area. 

B. STANDARD ERROR OF GPS OBSERVATIONS 

To determine the accuracy of GPS, measurements were taken at bench mark 
TREE on the NPS campus. Three measurements, each of about three hours 
duration, were taken for this analysis. Double difference fixed solutions were 
used to analyze the data. The double difference fixed solutions for these 
measurements are given in Table 30. 



Table 30. Ah AT STATION TREE 

r I X 



Date 


Time (UTC) 


Ahpix ( m ) 


Feb. 5, 88 


0730- 1025 


-8.4262 


Feb. 7, 88 


0721 - 1004 


-8.3532 


Feb. 11, 88 


0705 - 0959 


-8.3602 



The standard error of the mean of Ah (a / Vn) is about ± 2 cm at TREE and 
the standard error (a) of the single measurement is about ± 4 cm for three hours 
of observation. 

Results obtained by Strange in 1985 for a project southeast of Phoenix, 
Arizona, show that using multiple GPS observations at the same point lead to 
uncertainties typically less than 2 cm in Z-direction over 20-km distances 
[Zilkoski, 1988]. 

Thus the accuracy of the ellipsoid differences obtainable by GPS 
measurements is expected to be about ± 2 cm for about three hours of 
observation. 

C. CORRECTIONS FOR SURFACE METEOROLOGICAL 
VALUES 

The meteorological correction is a function of the atmospheric refractivity 
as computed by the surface meteorological values of pressure, temperature, 
humidity and the elevation angle of the satellite. Larger corrections are required 
for low elevation angles, as the signal travels a longer path through the 
troposphere. A modified Hopfield troposphere model [Fell, 1975] is used for 
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automatic processing of Trim640 software. The model corrects for at least 90% 
of the tropospheric delay [Remondi, 1984]. Default surface meteorological 
values for automatic processing are 1010 millibars, 20° C and 50 % relative 
humidity. These parameters can be reset before the processing. Trim640 allows 
only one pressure, temperature and humidity entry for each site per session. 
Ellipsoid height differences, Ahj^, were calculated for four marks on the NPS 

campus and two off campus using the default meteorological parameters and 
using the true values. There are given in Table 31 along with differences found 
between the two methods. 



Table 31. COMPARISON OF SURFACE METEOROLOGICAL 
CORRECTION 



Bench mark 


Default Ahp^ (m) 


True Ah FIX (m) 


Difference (mm) 


GH2 


-3.7777 


- 3.7766 


1.1 


GH3 


- 0.6376 


-0.6383 


0.6 


GH4 


0.9070 


0.9055 


1.5 


GH7 


-3.8524 


-3.8509 


1.5 


K 152 


2.2398 


2.2430 


-3.2 


J 697 


-71.3179 


-71.3023 


15.6 



The difference is about ± 2 mm on the NPS campus and about ± 2 cm in the 
Monterey Bay area. The Trim640 program contains no ionospheric model. For 
first-order relative positioning (1 part in 10 5 ) and a baseline shorter than 50 km 
ionospheric delay can be neglected [Remondi, 1984]. In this study surface 
meteorological corrections were made during data processing to get the best Ah 
solutions. 



D. OBSERVATION DURATIONS FOR GPS 

To find out how long observations must be made to obtain the best solutions 
in the field when batteries are used for power, the GPS data from S 812, J 697 
and K 152 were segmented into averaging periods of length nAt, where At = 10 
minutes and n = 1, 2, 3, ...,10. Default meteorological values were used for this 
processing. Ah FK for various observation durations for the stations is given in 

Table 32 and plotted in Figure 1 1 to Figure 13. 
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Table 32. Ah^ AS A FUNCTION OF OBSERVATION DURATION 

FIX 



Time (min) 


S 812 (m) 


J 697 (m) 


K 152 (m) 


10 


1.5309 


71.5291 


2.3135 


20 


1.6141 


70.9697 


2.3102 


30 


1.6441 


71.2185 


2.3062 


40 


1.6539 


71.4514 


2.3065 


50 


1.6558 


71.4762 


2.3060 


60 


1.6498 


71.5343 


2.2563 


70 


1.6464 


71.5309 


2.2504 


80 


1.6349 


71.3407 


2.2504 


90 




71.3179 


2.2548 


100 






2.2430 




Figure 11. Ah vs. observation durations at station S 812. 
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Figure 12. Ah vs. observation durations at station J 697. 




Figure 13. Ah vs. observation durations at station K 152. 

It is clear that Ah was affected by the observation duration. Ah varied quite 
a lot in the first forty minutes. After about forty minutes, Ah tended to close to 
the mean value of the observations. After sixty minutes, Ah didn't vary much. 
Trimble recommends, to ensure sufficient quality data and to obtain first-order 
results on single observations, at least one hour should be taken for a four or five 
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satellite session. For baseline lengths up to 30 km the conditions for the fixed 
double solution of automatic processing will be met if data are collected for 
ninety minutes on four or five satellites [Trimble, 1987b]. The results shows that 
the best solutions can be improved with longer observation times. Normally, at 
least sixty minutes of observations are required for vertical control. 



58 



VII. RESULTS AND DISCUSSION 



A. EVALUATATION OF GEOID MODEL 

To find the accuracy of the geoid model the orthometric heights of the 
check marks from GPS, Hq PS , were calculated, and these were compared to the 
orthometric height of the check marks from the levelling, H LEVELLING 

I started with the assumptions that C q = I and Cpq = 0. Five-parameter and 
six-parameter geoid models were used to calculate the orthometric height of the 
check marks. GWM 27 is a check mark in the off campus area. GH7 and GH8 
are the check marks on the NPS campus. N 0 (X,Y,Z) was calculated using both 

the NGS and the Trimvec programs. 

1. Hcp g from Five-Parameter Geoid Model 
a. Hqps in the off Campus Area 

H cps in GWM 27 was calculated by using seven control marks, six 
control marks (excluding B 21 or J 697) and five control marks (excluding B 21 
and J 697). Comparisons of H in GWM 27 are given in Table 33. 



Table 33. COMPARISONS OF H USING FIVE-PARAMETER 
MODEL AT GWM 27 



GWM 27 


HgPS ' ^LEVELLING (cm) 


# of control marks 


Trimvec 


NGS 


7 


-11.55 


-2.24 


6 (excluding B21) 


- 9.77 


-8.07 


6 (excluding J 697) 


-4.87 


-3.65 


5 


2.14 


6.58 



b. Hqps on the NPS Campus 

Hqpj in GH7 and GH8 were calculated by using seven control marks. 
Comparisons of H in GH7 and GH8 are given in Table 34. 
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Table 34. COMPARISONS OF H USING FIVE-PARAMETER 



MODEL AT GH7 AND GH8 



Five parameters 


^GPS ' ^LEVELLING ( Cm ) 


Check mark 


Trimvec 


NGS 


GH7 


0.52 




0.82 


GH8 


-0.20 




0.00 



2. H cpr from Six-Parameter Geoid Model 
a. Hqps in the off Campus Area 

H gps in GWM 27 were calculated by using seven control marks, and 
six control marks (excluding B 21 or J 697). Comparisons of H at GWM 27 are 
given in Table 35. 



Table 35. COMPARISONS OF H USING SIX-PARAMETER 
MODEL AT GWM 27 



GWM 27 


^cps _ ^levelling ( cm ) 


# of control marks 


Trimvec 


NGS 


7 


-12.40 


-11.13 


6 (excluding B21) 


- 9.37 


- 7.72 


6 (excluding J 697 


i -3.49 


-1.01 



b. Hqps on the NPS Campus 

Hqps of GH7 and GH8 were calculated by using seven control marks. 
Comparisons of H at GH7 and GH8 are given in Table 36. 



Table 36. COMPARISONS OF H USING SIX-PARAMETER 
MODEL AT GH7 AND GH8 



Six parameters 


^GPS " ^LEVELLING ( cm ) 


Check mark 


Trimvec 


NGS 


GH7 


1.22 




1.22 


GH8 


0.22 




0.25 



Since the accuracy of predicted geoid height at check marks (± 0.2 to 2 
cm) is less than the accuracy of GPS observations (± 2 to 4 cm), it may be 
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concluded that P = I is a satisfactory covariance matrix and further computation 
is unnecessary. 

B. ACCURACY 

The best accuracy of the five-parameter geoid model is ± 2 cm in the off 
campus area and is ± 2 to 10 mm for the NPS campus. The best accuracy of six- 
parameter geoid model is ± 1 cm in the off campus area and is ± 2 to 10 mm on 
the NPS campus. 

C. DISCUSSION 

The errors in the geoid model include errors associated with GPS and with 
levelling. 

1. Errors Associated With GPS 

a. Atmospheric Delays 

The time delay of RF signals passing through the ionosphere is due to 
a reduction in speed and the bending of the ray, both effects being due to 
refraction. The overall delay in the signal is nearly inversely proportional to the 
square of the frequency. The ionospheric delay may be calculated by using two- 
frequency receivers (C/A code, P code), and corrected for with a satisfactory 
degree of accuracy. Tropospheric delays are independent of frequency. They 
are relatively small and can be modelled using the elevation angle of the satellite. 
The Trimble 4000SX receiver only observes the C/A code . The Trim640 
program contains no ionospheric modelling. The modified Hopfield model [Fell, 
1975] is used for the tropospheric delay. The combined effects of unmodelled 
ionospheric and tropospheric errors are estimated to result in a satellite- to- user 
range error of from 2.4 to 5.2 m [Milliken, 1980]. 

b. Selection of Appropriate Control Marks 

Only eight permanent bench marks were recovered in the off campus 
study area. D 697 was obstructed by a tree. F 813 was set near a steel bridge. 
Although offset levelling was performed, the ellipsoid height could have been 
affected by offsetting. The shape of the bench marks were spread out roughly in 
a rhomboid shape. B 21 and J 697 are on the ends of the rhomboid. The 
distance is about 33 km from B 21 to J 697. The elevation changes form 5.787 to 
85.061 m. If more bench marks could be recovered and be connected to these 
marks, the accuracy of the seven control marks could be improved. 
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c. Error of Antenna Height Measurements 

The slant heights of GPS antenna were measured by using a steel tape 
before and after data collection. Thus, the heights of antenna were calculated. 
The systematic error of the tape, reading error, calculation error and a slack tape 
doing measurement will affect the accuracy of the ellipsoid height difference. 
The error of antenna height measurement is ± 2 to 10 mm. 

d. Selecting observation period with optimal satellite geometry 

Because there are only seven GPS satellites in operation at present, 

the observing window is only four to five hours a day for four or five satellites. 
The observation schedule should be planed early and considered the travel time, 
and time to reoccupy the marks. In this study most of the GPS measurements 
were made at night. PDOP were not always optimum. 

e. Error in the Computing of Trim640 Program 

There is no standard specification of the GPS software at present. 
Round off error in the data processing, the data record frequency, and the 
precision of the ephemeris could affect the components of the baseline. 
Specifications of the GPS software certainly will be developed in the future. 

2. Errors in Levelling 

The error in levelling would cause the errors of the observation 
equations when determining the local geoid model. The maximum allowable 
closing error between the forward and backward running of a section is 3.0 mm 
Vk for first-order class I levelling, and 9.0 mm v/k for third-order class I 
levelling, where k is the length of the section measured in km [Bodnar, 1975]. 
The average distance from K 152 to the permanent bench marks is about 13 km 
in the off campus study area. Thus, for first-order class I levelling, the 
maximum allowable closing error would be about 1 1 mm in the off campus area. 
The first-order bench marks were set around 1933. The elevations of the marks 
may have changed caused by subsidences or earthquakes. Since the time did not 
permit rerunning the levelling for the off campus area, published elevations were 
assumed correct. Errors in levelling could contribute to errors in the geoid 
model. 

The distance of the interlocking levelling circuit on the NPS campus is 
about 0.8 km. For third-order class I levelling, the maximum allowable closing 
error is about 7 mm on the NPS campus. The 5 mm closing error was obtained 
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for that region. The geoid height model thus had a better accuracy on the NPS 
campus. 
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VIII. CONCLUSIONS AND RECOMMENDATIONS 



Local geoid models for the Monterey, California, area were determined by 
GPS differential positioning using the method of collocation. The local geoid 
models were based on Rapp's 360 degree x 360 order global geoid models 
determined from gravity measurements. Control data were adjusted by least 
squares to solve for the parameters in the local geoid model. Also studied were 
factors that affected the GPS-measured ellipsoid height differences. These included 
(1) comparing GPS differencing solutions, (2) standard error of GPS observations, 
(3) corrections for surface meteorological values, and (4) observation durations for 
GPS. Conclusions are the following: 

1 . The accuracy of a local geoid height, N, is the same, whether the NGS or 
the Trimvec version of Rapp’s global geoid model is used, although N 0 (X,Y,Z) 

calculated using Trimble's Trimvec program gives the best results for points on the 
NPS campus, and N 0 (X,Y,Z) from the NGS program gives the best results for the 

larger study area. For practical use the orthometric heights can be calculated by 
utilizing the local geoid model. 

2. The areas where the geoid is smoothest lead to the best results. This was 
found by comparing the local geoid model used in the larger Monterey Bay area to 
the local geoid model used on the NPS campus. Also the density of control marks on 
the NPS campus was much higher than for the larger study area. If a higher 
accuracy is required in a large area where the geoid may be undulating, more 
control marks are required. These control marks should be strategically located 
throughout the network in order to determine the geoid slope. 

3 . The local geoid model does improve the accuracy of Rapp's global geoid 
model locally. The accuracy for Rapp’s model is ± 2 m in an area 55 km x 55 km. 
The accuracy for the local geoid model is ± 2 cm in an area 30 km x 30 km and ± 1 
cm in an area 0.5 km x 0.5 km. 

4. A GPS double difference fixed solution gives the best ellipsoid height 
differences, Ah, for baseline lengths up to 30 km. 

5. The standard error of GPS observations of the ellipsoid height difference, 
Ah, is ± 2 cm to ± 4 cm, so there is no need to perform the collocation with a non- 
unit weight matrix. Also it is not possible to reduce the error of the GPS 
observation below ± 2 cm. 
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6. There is no need for local surface meteorological corrections to the GPS 
observations. Default meteorological values (1010 millibars, 20° C and 50% 
relative humidity) are sufficient to meet the required accuracy. 

7. The longer the duration of the GPS measurements, the better the result will 
be. The accuracy of the geoid model for the NPS campus, where GPS 
measurements were made over four- to five-hour periods, is better than the 
accuracy of the geoid height model in the larger study area, where measurements 
were made for periods of only ninety minutes. 

8. A five-parameter geoid model gives the best solution for the NPS campus, 
while a six-parameter geoid model gives the best solution for the larger study area. 
The five-parameter geoid model which does not include a AZ term should be used 
with caution where AZ greatly exceeds 5 m. 

Recommendations are as follows: 

1 . The GPS antenna should be set at least 1 m above the ground to reduce the 
multiple effects from the ground. The height of the antenna should be carefully 
measured before and after the data collection. 

2. The configuration of the control marks should be carefully planned. Sites 
should be selected to optimize connections to stations with known precise 
orthometric height differences, minimizing the effects of obstructions. Control 
marks in a large area should be numerous enough to meet the required accuracy. 

3. To minimize computation errors least squares adjustments should be made 
using double precision processing when solving the geoid model. For large 
amounts of GPS data a mainframe computer should be used rather than a 
microcomputer to reduce processing time. 
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APPENDIX A. FORTRAN PROGRAM GPSCON 



* WRITTEN BY MA, WEI-MING 05/08/88 

* THIS PROGRAM READS DATA FROM GEOID DATA FILE AND CALLS 

* THE SUBROUTINE CONTOR TO PLOT GEOID HEIGHTS. 

* THE VARIABLES USED ARE: 

* GEOHT : GEOID HEIGHTS (m) 

* NX : NUMBER OF POINTS IN THE X-DIRECTION 

* NY : NUMBER OF POINTS IN THE Y-DIRECTION 

* 

REALM GEOHT(22,15) 

DATA NX,NY/22,15/ 

C 

C READ GEOID HEIGHTS FROM GEOID 1 DATA FILE 
C 

CALL EXCMS ('FILED EF 8 DISK GEOID1 DATA Al’) 

DO 200 J = 1, NY 

DO 1001 = NX, 1,-1 

READ(8,*) GEOHT (I, J) 

100 CONTINUE 
200 CONTINUE 
C 

C PLOT THE GEOID HEIGHT BY CALLING SUBROUTINE CONTOR 
C 

CALL CONTOR(GEOHT,NX,NY) 

STOP 

END 

SUBROUNTINE CONTOR(A,NX,NY) 

* * 

* WRITTEN BY MA, WEI-MING 05/08/88 

* THIS SUBROUTINE CONTOURS AN NX BY NY ARRAY OF 

* REGULARLY SPACED POINTS. 

* NOTE: THIS ARRAY MUST BE REALM. 

* THE VARIABLES USED ARE: 

* A : SINGLE PRECISION NX BY NY ARRAY OF 

* REGULARLYSPACED POINTS. 

* NX : NUMBER OF POINTS IN THE X-DIRECTION 

* NY : NUMBER OF POINTS IN THE Y-DIRECTION 

* ZINC : CONTOUR INTERVAL 

* 

****************************%********%**%**********%***************%** 

DIMENSION A(NX,NY) 

COMMON WORK(5000) 

C 

C SET PARAMETERS FOR AXES 
C 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 



★ 

* 

* 

★ 

* 

* 

* 

* 
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no non non non non 



XORIG = -55.0 
XSTP = 2.0 
XMAX = -34.0 
YORIG = 35.0 
YSTP = 2.0 
YMAX = 49.0 

SET CONTOUR LEVEL 



ZINC = 0.1 
CALL COMPRS 
CALL SETCLR('CYAN') 

SET PAGE AND PLOT SIZES, SET UP AXES FOR PLOT 

CALL PAGE(8. 0,6.0) 

CALL BCOMON(5000) 

CALL AREA2D(7 .0,5.0) 

LABEL AXES 

CALL XNAME('LONGITUDE 121 DEGREE WEST (MINUTES)$', 100) 
CALL YNAME('LATTTUDE 36 DEGREE NORTH (MINUTES )$’, 100) 

CALL GRAF(XORIG,XSTP,XMAX, Y ORIG,Y STP, YMAX) 

CALL FRAME 

MAKE CONTOURS AND DRAW 

CALL SETCLR(’RED’) 

CALL CONMIN(3.0) 

CALL CONANG(60.) 

CALL CONLIN(0,'MYCON7LABELS', 1,10) 

CALL CONMAK(A,NX,NY,ZINC) 

CALL CONTUR(l, 'LABELS', DRAW') 

END PLOT 

CALL ENDPL(O) 

CALL DONEPL 
RETURN 
END 

C SUBROUTINE MY CON(RARA Y,IARAY) 

********************************************************************** 

* * 

* THIS SUBROUTINE MAKES NEGATIVE CONTOURS DASHED AND * 

* THE ZERO LINE HEAVYIER * 

* * 

********************************************************************** 

DIMENSION RARA Y (2),I ARA Y ( 1 ) 

CALL RESET('DASH') 

IF (RARAY(l).GE. 0.) GO TO 10 
CALL DASH 
10 RARAY(2) = 1. 

IARAY(l) =1 
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IF (RARAY(l).EQ.O.) IARAY(1) = 2 

RETURN 

END 
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APPENDIX B. BASIC PROGRAM LOBS.BASIC 



5 REM LEAST SQUARES BY OBSERVATION 
10 DIM L(N%,1) 

20 DIM A(N%,X%), P(X%,X%), EX(X%,1) 

30 DIM CX(X%,1) 

40 DIM AT(X%,N%), R(50), ATP(X%,X%) 

50 DIM ATPA(X%,X%), ATPL(X%,1), AI(X%,X%) 

69 PRINT "INPUT# OF OBSERVATIONS" 

70 NPUT NO 

200 PRINT "INPUT# OF PARAMETERS" 

270 INPUT N1 

350 ERASE L 

351 ERASE A 

352 ERASE AT 

353 ERASE P 

354 ERASE ATP 

355 ERASE ATPL 

356 ERASE AI 

357 ERASE CX 

358 ERASE EX 

359 ERASE ATP A 

380 DIM L(NO,l), A(NO,Nl), AT(Nl,NO), P(NO,NO), ATP(Nl,NO), 

ATPL(N 1 , 1 ), AI(N 1 ,N 1 ) 

390 DIM CX(N1,1), EX(NO,l), ATPA(Nl.Nl) 

400 ITERO = 0 

1 1 80 PRINT " INPUT COEFFICIENTS OF OBSERVATION EQUATIONS AND 
WEIGHTS" 

1190 FOR K = 1 TO NO 
1200 FOR J = 1 TO N1 
1220 PRINT "COEFFICIENT ", K, J 
1230 INPUT A(K,J) 

1240 PRINT A(K,J) 

1250 NEXT J 

1 260 PRINT " INPUT OBSERVED VALUE ", K 
1270 INPUT L(K,1) 

1280 PRINT " INPUT WEIGHT OF OBSERVED VALUE ", K 
1300 INPUT P(K,K) 

1390 NEXT K 
1410 DOF = 0 

1420 REM LEAST SQUARES 
1430 M = NO 
1450 L = N1 

1460 PRINT " M =", M, "N1 = ", Nl, "L = "L 
1470 FOR 1 = 1 TO M 
1480 FOR J = 1 TO L 
1490 AT(J,I) = A(I,J) 

1500 NEXT J 

1510 NEXT I 

1520 REM ATP = AT * P 

1530 N = NO 

1540 FOR 1=1 TO L 
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1550 FOR J =1 TON 

1560 ATP(I,J) = 0 

1570 FOR K = 1 TO M 

1580 ATP(I,J) = ATP(I,J) + AT(I,J) * P(K,J) 

1590 NEXT K 
1600 NEXT J 
1610 NEXT I 

1620 REM ATPA = ATP * A 
1630 N = N1 
1640 FOR I = 1 TO L 
1650 FOR J = 1 TO N 
1660 ATPA(I,J) = 0 
1680 ATPA(IJ) = 0 
1690 FOR K = 1 TO M 

1710 ATPA(I,J) = ATPA(1,J) + ATP(I,K) * A(K,J) 
1720 NEXT K 
1730 NEXT J 
1740 NEXT I 

1750 REM ATPL = ATP * L 
1760 N = 1 
1770 FOR I = 1 TO L 
1780 FOR J = 1 TO N 
1790 ATPL(I,J) = 0 
1800 FOR K = 1 TOM 

1810 ATPL(I,J) = ATPL(I,J) + ATP(I,K) * L(K,J) 
1820 NEXT K 
1830 NEXT J 
1840 NEXT I 

1970 REM AI = INV(ATPA) 

1980 I = N1 

1990 M = N1 

2000 N = I - 1 

2020 MI = M - 1 

2030 FOR J = 1 TO I 

2040 FOR K = 1 TO I 

2050 AI(J,K) = ATPA(J,K) 

2055 NEXT K 

2060 NEXT J 

2070 FOR K =1 TO I 

2080 FOR J = 1 TO MI 

2090 R(J) = AI(1,J+1) / AI(1,1) 

2100 NEXT J 

2110 R(M) = 1! / AI(1,1) 

2120 FOR L = 1 TO N 
2130 FOR J = 1 TO MI 

2140 AI(I,J) = AI(L+1,J+1) - AI(L+1,1) * R(J) 
2150 NEXT J 

2160 AI(L,M) = -AI(L+1,1) * R(M) 

2170 NEXT L 
2180 FOR J = 1 TO M 
2190 AI(I,J) = R(J) 

2200 NEXT I 
2210 NEXT K 
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2220 

2230 

2240 

2250 

2260 

2270 

2280 

2300 

2310 

2320 

2325 

2330 

2340 

2480 

2490 

2500 

2510 

2520 

2530 

2540 

2560 

2570 

2580 

2590 

2600 

2610 

2620 

2630 

2640 

2650 

2660 

2670 

2680 

2690 

2695 

2700 

2710 

2720 



REM CX = AI * ATPL 
L = N1 
N = 1 
M = N1 

FOR 1=1 TOL 
FOR J = 1 TO N 
CX(I,J) = O 
FOR K = 1 TO M 

CX(I,J) = CX(I,J) + AI(I,K) * ATPL (K,J) 
NEXT K 

PRINT "PARAMETER # ", I, CX(I,J) 
NEXT J 
NEXT I 

PRINT " RESIDUAL" 

L = NO 
N = 1 
M = N1 

REM EX = A * CX 
FOR I = 1 TO L 
FOR J = 1 TO N 
EX(I,J) = 0 
FOR K = 1 TO M 

EX(I,J) = EX(I,J) + A(I,K) * CX(K,J) 

NEXT K 

NEXT J 

NEXT I 

SD = 0 

FOR K = 1 TO NO 
V = EX(K,1) - L(K,1) 

PRINT K, V 
SD = SD + V * V 
NEXT K 

SD = SQR(SD / (NO - N1 + DOF)) 

PRINT SD 

PRINT " STD.DEV", SD 
ITER = ITER + 1 
IF ITER < 3 GOTO 350 
END 
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APPENDIX C. FORTRAN PROGRAM GPSDIS 






* 

* 

* 

* 

* 



WRITTEN BY WEI-MING MA 05/14/88 

THIS PROGRAM READS DATA FROM TERMINAL AND COMPUTES 
THE ORTHOMETRIC HEIGHTS FOR THE CHECK MARKS THEN 
COMPUTES THE DIFFERENCE. 

THE VARIABLES USED ARE: 



* 


RX,RY,RZ 


THE 


* 


CX,CY,CZ 


THE 


* 


DX,DY,DZ 


THE 


* 


H 


THE 


* 


LH 


THE 


* 


DH 


THE 


* 


GN 


THE 


* 


LH 


THE 


* 


DIF 


THE 


* 


A,B,C,D 


THE 


* 


A1,B1,C1 


THE 


* 


D1,E1 




* 


IPA, IPA 1 : 


THE 


* 


IW 


THE 



* 

* 

* 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 



REAL* 8 RX, RY, RZ, CX, CY, CZ, DX, DY, DZ, LH, GN, 

REAL*8 A, B, C, D, Al, Bl, Cl, Dl, El, H, DH, DIF, 

INTEGER IPA, IPA1.IW 

CHARACTER* 1 RESPN1, RESPN2, RESPN3 

IW = 6 

IPA = 0 

IPA1 =0 



READ DATA FROM TERMINAL 



10 PRINT *, 'ENTER # OF PARAMETERS' 

READ * IPA 

IF (IPA ^EQ. 5 .OR. IPA .EQ. 6) THEN 

PRINT *, 'ENTER THE X COORDINATE IN WGS 84 OF THE ', 
& 'REFERENCE MARK' 

READ *, RX 

PRINT *, 'ENTER THE Y COORDINATE IN WGS 84 OF THE ’, 
& 'REFERENCE MARK' 

READ *, RY 

PRINT *, 'ENTER THE Z COORDINATE IN WGS 84 OF THE ', 
& 'REFERENCE MARK' 

READ *, RZ 
ELSE 
STOP 
END IF 
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IF (RESPN1 .EQ. 'Y') THEN 

PRINT *,'D0 YOU WANT TO CHANGE THE CHECH MARK ("Y" OR 
& '"N”)' 

READ *, RESPN3 
IF (RESPN3 .EQ. 'Y') GO TO 30 
GO TO 40 
END IF 

20 PRINT *,'DO YOU WANT TO CHANGE THE REFERENCE MARKS ("Y” ', 
& 'OR "N")’ 

READ *, RESPN1 

IF (RESPN1 .EQ. 'Y') GO TO 10 

30 PRINT *, 'ENTER THE X COORDINATE IN WGS 84 OF THE CHECK.’, 

& 'MARK' 

READ * CX 

PRINT *,' 'ENTER THE Y COORDINATE IN WGS 84 OF THE CHECK MARK' 
READ * CY 

PRINT VENTER THE Z COORDINATE IN WGS 84 OF THE CHECK MARK' 
READ *, CZ 

PRINT *, 'ENTER THE ELLIPSOID DIFFERENCE, DH’ 

READ *, DH 

PRINT *, 'ENTER THE GLOBAL GEOID HEIGHT OF THE CHECK MARK’, 

& ', GN' 

READ *, GN 

PRINT *, 'ENTER THE LEVELED OTHOMETRIC HEIGHT OF THE CHECK ’ 
& , 'MARK' 

READ *, LH 
DX = CX - RX 
DY = CY - RY 
DZ = CZ- RZ 

IF (IPA .EQ. IPA1 .AND. RESPN1 .EQ. ’N' .AND. IPA .EQ. 5) GO TO 50 
IF (IPA .EQ. EPA1 .AND. RESPN1 .EQ. 'N' .AND. IPA .EQ. 6) GO TO 60 
40 IF (EPA .EQ. 5) THEN 

PRINT *, 'ENTER THE ELLIPSOID HEIGHT HO' 

READ *, HO 

PRINT *, 'ENTER THE COEFFICIENT A’ 

READ *, A 

PRINT *, 'ENTER THE COEFFICIENT B' 

READ *, B 

PRINT *, 'ENTER THE COEFFICIENT C' 

READ * C 

PRINT *, 'ENTER THE COEFFICIENT D' 

READ *, D 

THE 5-PARAMETER GEOID MODEL 

50 H = -GN + DH + HO + A*DY + B*DX**2 + C*DY**2 + D*DX*DY 

ELSE 

PRINT *, 'ENTER THE ELLIPSOID HEIGHT HO’ 

READ *, HO 

PRINT *, 'ENTER THE COEFFICIENT A'" 

READ * A1 

PRINT *, 'ENTER THE COEFFICIENT B'" 

READ *, B1 
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PRINT *, 'ENTER THE COEFFICIENT C m 
READ * Cl 

PRINT *, 'ENTER THE COEFFICIENT D’" 

READ *, D1 

PRINT *, 'ENTER THE COEFFICIENT E"' 

READ *, El 

THE 6-PARAMETER GEOID MODEL 

60 H = -GN + DH + HO + A1*DX + B1*DZ + C1*DX**2 + D1*DY**2 
& + E1*DX*DY 

END IF 

COMPUTE THE DIFFERENCE AND PRINT THE RESULT 

DIF = ( H - LH) * 100. 

WRITE(IW, 1) IPA 
WRITE(IW, 2) LH, H, DIF 

1 FORMAT (5X, II,' PARAMETERS') 

2 FORMAT (/5X/H-LEVELING =’, F10.3, ' M' 

& /5X/H-GPS =', F10.3, ' M' 

& /5X, ’DIFFERENCE =’, F10.3, ' CM’/) 

IPA 1 = IPA 
H = 0. 

DIF = 0. 

PRINT *, 'MORE COMPUTATION ("Y" OR "N”)’ 

READ *, RESPN2 

IF (RESPN2 .EQ. ’Y' ) GO TO 20 

STOP 

END 
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APPENDIX D. FORTRAN PROGRAM DISTCO 






* WRITTEN BY MA, WEI_MING 05/11/88 * 

* THIS PROGRAM COMPUTES THE DISTANCE BETWEEN THE BENCH * 

* MARKS THEN CALLS THE SUBROUTINE C TO CALCULATE THE * 

* COVARIANCE BETWEEN THE BENCH MARKS. * 

* THE VARIABLES USED ARE: * 

* X,Y,Z : THE COORDINATES OF X,Y,Z, IN WGS 84 (m) * 

* DIS : THE DISTANCES BETWEEN THE MARKS (m) * 

* CQ : THE COVARIANCES BETWEEN THE BENCH MARKS * 



* * 
********************************************************************** 



REAL*8 X(10), Y(10), Z(10) 

REAL* 8 DIS(10,10), CQ(10,10) 

DATA DIS/100*0./, CQ/100*0./ 

IW = 9 
N = 7 
C 

C READ X, Y, Z FROM THE DATA FILE 
C 

CALL EXCM S (’FILEDEF 8 DISK MBXYZ DATA Al') 

DO 1001= 1,N 

READ(8,*) X(I), Y(I), Z(I)) 

100 CONTINUE 
C 

C COMPUTE THE DISTANCE AND PRINT THE RESULTS 
C 

DO 300 I = 1, N 
WRITE(IW,1) I 

1 FORMAT(/lX,’FROM BENCH MARK #’, 12,' TO’, 'DISTANCE (m)’/ 

& IX, 53(’-')) 

DO 200 J = 1, N 

IF (J .EQ. I) GO TO 10 

DIS(I,J) = SQRT((X(J) - X(I)) ** 2 + (Y(J) - Y(I)) ** 2 
& +(Z(J)-Z(I))**2) 

10 WRITE(IW,2) J, DIS(I,J) 

2 FORMAT(l IX, 12, 14X, G20.8) 

200 CONTINUE 

300 CONTINUE 
C 

C CALL SUBROUTINE C TO CALCULATE THE CQ AND PRINT THE RESULTS 
C 

WRITE (TW,3) 

3 FORMAT(/3X, ’CQ = ’) 

CALL C(DIS,N,CQ) 

DO 400 J = 1, N 

WRITE(IW,4) (CQ(I,J), I = 1, N) 

4 FORMAT(3X, 7(G16.7, 2X) 

400 CONTINUE 
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STOP 

END 

^SUBROUNTINE C(DIS,N,CQ) 

* * 

WRITTEN BY MA, WEI_MING 05/11/88 * 

THIS SUBROUTINE COMPUTES THE COVARIANCES OF THE RANDOM * 
ERROR DISTRIBUTION FUNCTION: * 

C(R) = A + B * SIN(D*DIS) * 

THE VARIABLES USED ARE: * 

A : STANDARD DEVIATION OF THE OBSERVATIONS * 

B : CONSTANT * 

C : CONSTANT (RADIANS/METERS) * 

DIS : DISTANCE (METERS) * 

* 

REAL*8 DIS(10,10), CQ(IO.IO), A, B, D 
A = 0.137153 
B = 0.147951 
D = 2.85958 
DO 200 J = 1 , N 
DO 1001= 1,N 

CQ(I,J) = A - B * SIN(D * DIS(I,J) / 33300.) 

100 CONTINUE 
200 CONTINUE 
RETURN 
END 



* 

* 

* 

* 

* 

* 

* 

* 

* 

★ 
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